Changeset 22898

Timestamp:

07/04/18 07:48:56 (7 years ago)

Author:

bdef

Message:

NEW:deactivation of hydrologydc for frozen elements

Location:

issm/trunk-jpl/src/c

Files:

• analyses/HydrologyDCInefficientAnalysis.cpp (modified) (7 diffs)
• analyses/HydrologyDCInefficientAnalysis.h (modified) (3 diffs)
• classes/Elements/Element.cpp (modified) (1 diff)
• classes/FemModel.cpp (modified) (169 diffs)
• classes/FemModel.h (modified) (5 diffs)
• cores/hydrology_core.cpp (modified) (4 diffs)
• cores/transient_core.cpp (modified) (11 diffs)
• shared/Enum/EnumDefinitions.h (modified) (1 diff)
• shared/Enum/EnumToStringx.cpp (modified) (1 diff)
• shared/Enum/StringToEnumx.cpp (modified) (7 diffs)
• solutionsequences/solutionsequence_hydro_nonlinear.cpp (modified) (12 diffs)

issm/trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.cpp

@@ -100 +100 @@
         iomodel->FetchDataToInput(elements,"md.initialization.sediment_head",SedimentHeadHydrostepEnum);
         iomodel->FetchDataToInput(elements,"md.hydrology.sediment_transmitivity",HydrologydcSedimentTransmitivityEnum);
+        iomodel->FetchDataToInput(elements,"md.hydrology.mask_thawed_node",HydrologydcMaskThawedNodeEnum);
         if(iomodel->domaintype!=Domain2DhorizontalEnum){
                 iomodel->FetchDataToInput(elements,"md.mesh.vertexonbase",MeshVertexonbaseEnum);
@@ -184 +185 @@
 
         /*Intermediaries*/
+        bool     thawed_element;
         int      domaintype;
         Element* basalelement;
@@ -199 +201 @@
                 default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
         }
+
+        Input* thawed_element_input = basalelement->GetInput(HydrologydcMaskThawedEltEnum); _assert_(thawed_element_input);
+        thawed_element_input->GetInputValue(&thawed_element);
+
+        /*Check that all nodes are active, else return empty matrix*/
+        if(!thawed_element) {
+        if(domaintype!=Domain2DhorizontalEnum){
+                        basalelement->DeleteMaterials();
+                        delete basalelement;
+                }
+                return NULL;
+        }
+
 
         /*Intermediaries */
@@ -297 +312 @@
 
         /*Intermediaries*/
-        int      domaintype;
+        bool             thawed_element;
+        int                      domaintype;
         Element* basalelement;
 
@@ -311 +327 @@
                         break;
                 default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
+        }
+
+        Input* thawed_element_input = basalelement->GetInput(HydrologydcMaskThawedEltEnum); _assert_(thawed_element_input);
+        thawed_element_input->GetInputValue(&thawed_element);
+
+        /*Check that all nodes are active, else return empty matrix*/
+        if(!thawed_element) {
+        if(domaintype!=Domain2DhorizontalEnum){
+                        basalelement->DeleteMaterials();
+                        delete basalelement;
+                }
+                return NULL;
         }
 
@@ -626 +654 @@
                 _error_("UnconfinedFlag is 0 or 1");
         }
-
-        /*Let's deal with the frozen parts for which transmitivity is zero*/
-        if (isthermal){
-                element->GetInputValue(&meltingrate,gauss,BasalforcingsGroundediceMeltingRateEnum);
-                element->GetInputValue(&groundedice,gauss,MaskGroundediceLevelsetEnum); 
-                if ((meltingrate<=0.0) && (groundedice>0)) sediment_transmitivity=0.;
-        }
-
         return sediment_transmitivity;
 }/*}}}*/
@@ -738 +758 @@
         }
 }/*}}}*/
+
+void  HydrologyDCInefficientAnalysis::HydrologyIDSGetMask(Vector<IssmDouble>* vec_mask, Element* element){
+        bool        active_element;
+        int         domaintype;
+        Element*    basalelement=NULL;
+
+        /*Get basal element*/
+        element->FindParam(&domaintype,DomainTypeEnum);
+        switch(domaintype){
+                case Domain2DhorizontalEnum:
+                        basalelement = element;
+                        break;
+                case Domain3DEnum:
+                        if(!element->IsOnBase()) return;
+                        basalelement = element->SpawnBasalElement();
+                        break;
+                default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
+        }
+        /*Intermediaries*/
+        int                                             numnodes    =   basalelement->GetNumberOfNodes();
+        IssmDouble*             meltingrate =   xNew<IssmDouble>(numnodes);
+        IssmDouble*             groundedice =   xNew<IssmDouble>(numnodes);
+
+        basalelement->GetInputListOnVertices(&meltingrate[0],BasalforcingsGroundediceMeltingRateEnum);
+        basalelement->GetInputListOnVertices(&groundedice[0],MaskGroundediceLevelsetEnum);
+
+        /*if melting rate is not positive and node is not floating, deactivate*/
+        for(int i=0;i<numnodes;i++){
+                if ((meltingrate[i]<=0.0) && (groundedice[i]>0)){
+                        vec_mask->SetValue(basalelement->nodes[i]->Sid(),0.,INS_VAL);
+                }
+                else{
+                        vec_mask->SetValue(basalelement->nodes[i]->Sid(),1.,INS_VAL);
+                }
+        }
+
+        if(domaintype!=Domain2DhorizontalEnum){
+                basalelement->DeleteMaterials();
+                delete basalelement;
+        }
+        xDelete<IssmDouble>(meltingrate);
+        xDelete<IssmDouble>(groundedice);
+}
+void HydrologyDCInefficientAnalysis::ElementizeIdsMask(FemModel* femmodel){/*{{{*/
+
+        bool     element_active;
+        Element* element=NULL;
+
+        for(int i=0;i<femmodel->elements->Size();i++){
+                element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
+                        Input* node_mask_input = element->GetInput(HydrologydcMaskThawedNodeEnum); _assert_(node_mask_input);
+
+                if(node_mask_input->Max()>0.){
+                        element_active = true;
+                }
+                else{
+                        element_active = false;
+                }
+                element->AddInput(new BoolInput(HydrologydcMaskThawedEltEnum,element_active));
+        }
+}/*}}}*/

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

@@ -1 @@
-/*! \file HydrologyDCInefficientAnalysis.h 
+/*! \file HydrologyDCInefficientAnalysis.h
  *  \brief: header file for generic external result object
  */
@@ -8 +8 @@
 /*Headers*/
 #include "./Analysis.h"
-class Node; 
+class Node;
 class Input;
 class HydrologyDCInefficientAnalysis: public Analysis{
@@ -40 +40 @@
                 IssmDouble GetHydrologyPVectorTransfer(Element* element, Gauss* gauss, Input* epl_head_input);
                 void ElementizeEplMask(FemModel* femmodel);
+                void HydrologyIDSGetMask(Vector<IssmDouble>* vec_mask, Element* element);
+                void ElementizeIdsMask(FemModel* femmodel);
 };
 #endif

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

@@ -2069 +2069 @@
                                 name==HydrologydcEplThicknessHydrostepEnum ||
                                 name==HydrologydcMaskEplactiveNodeEnum ||
+                                name==HydrologydcMaskThawedNodeEnum ||
                                 name==HydrologyHeadEnum ||
                                 name==HydrologyHeadOldEnum ||

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

@@ -136 +136 @@
         bool traceon=true;
         this->profiler=NULL; /*avoid leak, as we are not using the profiler ever in ad control run. */
-        
-        /*Store the communicator, but do not set it as a global variable, as this has already 
+
+        /*Store the communicator, but do not set it as a global variable, as this has already
          * been done by the FemModel that called this copy constructor: */
         IssmComm::SetComm(incomm);
@@ -152 +152 @@
         this->amrbamg = NULL;
         #endif
-        
+
         /*Save communicator in the parameters dataset: */
         this->parameters->AddObject(new GenericParam<ISSM_MPI_Comm>(incomm,FemModelCommEnum));
@@ -182 +182 @@
         if(parameters)delete parameters;
         if(results)delete results;
-        
+
         #if defined(_HAVE_NEOPZ_) && !defined(_HAVE_ADOLC_)
         if(amr)delete amr;
@@ -206 +206 @@
         /*First, recover the name of the restart file: */
         parameters->FindParam(&restartfilename,RestartFileNameEnum);
-        
+
         /*Open file for writing: */
         restartfid=pfopen(restartfilename,"wb");
@@ -215 +215 @@
         /*Create buffer to hold marshalled femmodel: */
         this->Marshall(NULL,&femmodel_size,MARSHALLING_SIZE);
-        femmodel_buffer=xNew<char>(femmodel_size); 
+        femmodel_buffer=xNew<char>(femmodel_size);
         /*Keep track of initial position of femmodel_buffer: */
         femmodel_buffer_ini=femmodel_buffer;
-        
+
         /*Marshall:*/
         this->Marshall(&femmodel_buffer,NULL,MARSHALLING_FORWARD);
@@ -376 +376 @@
 }/*}}}*/
 void FemModel::InitFromFids(char* rootpath, FILE* IOMODEL, FILE* toolkitsoptionsfid, int in_solution_type, bool trace, IssmPDouble* X){/*{{{*/
-        
+
         /*Initialize internal data: */
         this->solution_type    = in_solution_type;
         this->analysis_counter = nummodels-1;   //point to last analysis_type carried out.
         this->results          = new Results(); //not initialized by CreateDataSets
-        
+
         /*create IoModel */
         IoModel* iomodel = new IoModel(IOMODEL,in_solution_type,trace,X);
@@ -399 +399 @@
                 if(VerboseMProcessor()) _printf0_("      configuring element and loads\n");
                 ConfigureObjectsx(elements, loads, nodes, vertices, materials,parameters);
-                
+
                 if(i==0){
                         if(VerboseMProcessor()) _printf0_("      creating vertex PIDs\n");
-                        VerticesDofx(vertices,parameters); 
+                        VerticesDofx(vertices,parameters);
 
                         if(VerboseMProcessor()) _printf0_("      detecting active vertices\n");
@@ -409 +409 @@
 
                 if(VerboseMProcessor()) _printf0_("      resolving node constraints\n");
-                SpcNodesx(nodes,constraints,parameters,analysis_type_list[i]); 
+                SpcNodesx(nodes,constraints,parameters,analysis_type_list[i]);
 
                 if(VerboseMProcessor()) _printf0_("      creating nodal degrees of freedom\n");
@@ -487 +487 @@
         FILE* restartfid=NULL;
         char* restartfilename = NULL;
-        int   femmodel_size=0; 
-        int   fread_return=0; 
+        int   femmodel_size=0;
+        int   fread_return=0;
         char* femmodel_buffer=NULL;
         char* femmodel_buffer_ini=NULL;
@@ -513 +513 @@
 
         /*Figure out size of buffer to be read: */
-        fseek(restartfid, 0L, SEEK_END); 
+        fseek(restartfid, 0L, SEEK_END);
         femmodel_size = ftell(restartfid);
         fseek(restartfid, 0L, SEEK_SET);
 
         /*Allocate buffer: */
-        femmodel_buffer=xNew<char>(femmodel_size); 
+        femmodel_buffer=xNew<char>(femmodel_size);
 
         /*Read buffer from file: */
@@ -539 +539 @@
 void FemModel::SetCurrentConfiguration(int configuration_type,int analysis_type){/*{{{*/
 
-        /*Use configuration_type to setup the analysis counter, the configurations of objects etc ... but use 
-         * analysis_type to drive the element numerics. This allows for use of 1 configuration_type for several 
-         * analyses. For example: do a SurfaceSlopeX, SurfaceSlopeY, BedSlopeX and BedSlopeY analysis using the 
+        /*Use configuration_type to setup the analysis counter, the configurations of objects etc ... but use
+         * analysis_type to drive the element numerics. This allows for use of 1 configuration_type for several
+         * analyses. For example: do a SurfaceSlopeX, SurfaceSlopeY, BedSlopeX and BedSlopeY analysis using the
          * Slope configuration.*/
         int found=-1;
@@ -700 +700 @@
                         analyses_temp[numanalyses++]=EsaAnalysisEnum;
                         break;
-                
+
                 case SealevelriseSolutionEnum:
                         analyses_temp[numanalyses++]=SealevelriseAnalysisEnum;
@@ -825 +825 @@
 
         /*run solution core: */
-        profiler->Start(CORE);   
-        solutioncore(this); 
+        profiler->Start(CORE);
+        solutioncore(this);
         profiler->Stop(CORE);
 
         /*run AD core if needed: */
-        profiler->Start(ADCORE); 
-        ad_core(this);      
+        profiler->Start(ADCORE);
+        ad_core(this);
         profiler->Stop(ADCORE);
 
@@ -1121 +1121 @@
         /*Broadcast and plug into response: */
         ISSM_MPI_Allreduce ( &cpu_found,&cpu_found,1,ISSM_MPI_INT,ISSM_MPI_MAX,IssmComm::GetComm());
-        ISSM_MPI_Bcast(&response,1,ISSM_MPI_DOUBLE,cpu_found,IssmComm::GetComm()); 
+        ISSM_MPI_Bcast(&response,1,ISSM_MPI_DOUBLE,cpu_found,IssmComm::GetComm());
 
         /*Assign output pointers: */
@@ -1268 +1268 @@
         num_segments = mdims_array[counter-1];
 
-        /*Go through segments, and then elements, and figure out which elements belong to a segment. 
+        /*Go through segments, and then elements, and figure out which elements belong to a segment.
          * When we find one, use the element to compute the mass flux on the segment: */
         for(i=0;i<num_segments;i++){
@@ -1315 +1315 @@
         /*Figure out maximum across the cluster: */
         ISSM_MPI_Reduce(&maxabsvx,&node_maxabsvx,1,ISSM_MPI_DOUBLE,ISSM_MPI_MAX,0,IssmComm::GetComm() );
-        ISSM_MPI_Bcast(&node_maxabsvx,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());   
+        ISSM_MPI_Bcast(&node_maxabsvx,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
         maxabsvx=node_maxabsvx;
 
@@ -1339 +1339 @@
         /*Figure out maximum across the cluster: */
         ISSM_MPI_Reduce(&maxabsvy,&node_maxabsvy,1,ISSM_MPI_DOUBLE,ISSM_MPI_MAX,0,IssmComm::GetComm() );
-        ISSM_MPI_Bcast(&node_maxabsvy,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());   
+        ISSM_MPI_Bcast(&node_maxabsvy,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
         maxabsvy=node_maxabsvy;
 
@@ -1363 +1363 @@
         /*Figure out maximum across the cluster: */
         ISSM_MPI_Reduce(&maxabsvz,&node_maxabsvz,1,ISSM_MPI_DOUBLE,ISSM_MPI_MAX,0,IssmComm::GetComm() );
-        ISSM_MPI_Bcast(&node_maxabsvz,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());   
+        ISSM_MPI_Bcast(&node_maxabsvz,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
         maxabsvz=node_maxabsvz;
 
@@ -1406 +1406 @@
         /*Figure out maximum across the cluster: */
         ISSM_MPI_Reduce(&maxvel,&node_maxvel,1,ISSM_MPI_DOUBLE,ISSM_MPI_MAX,0,IssmComm::GetComm() );
-        ISSM_MPI_Bcast(&node_maxvel,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());   
+        ISSM_MPI_Bcast(&node_maxvel,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
         maxvel=node_maxvel;
 
@@ -1430 +1430 @@
         /*Figure out maximum across the cluster: */
         ISSM_MPI_Reduce(&maxvx,&node_maxvx,1,ISSM_MPI_DOUBLE,ISSM_MPI_MAX,0,IssmComm::GetComm() );
-        ISSM_MPI_Bcast(&node_maxvx,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());   
+        ISSM_MPI_Bcast(&node_maxvx,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
         maxvx=node_maxvx;
 
@@ -1454 +1454 @@
         /*Figure out maximum across the cluster: */
         ISSM_MPI_Reduce(&maxvy,&node_maxvy,1,ISSM_MPI_DOUBLE,ISSM_MPI_MAX,0,IssmComm::GetComm() );
-        ISSM_MPI_Bcast(&node_maxvy,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());   
+        ISSM_MPI_Bcast(&node_maxvy,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
         maxvy=node_maxvy;
 
@@ -1478 +1478 @@
         /*Figure out maximum across the cluster: */
         ISSM_MPI_Reduce(&maxvz,&node_maxvz,1,ISSM_MPI_DOUBLE,ISSM_MPI_MAX,0,IssmComm::GetComm() );
-        ISSM_MPI_Bcast(&node_maxvz,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());   
+        ISSM_MPI_Bcast(&node_maxvz,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
         maxvz=node_maxvz;
 
@@ -1502 +1502 @@
         /*Figure out minimum across the cluster: */
         ISSM_MPI_Reduce(&minvel,&node_minvel,1,ISSM_MPI_DOUBLE,ISSM_MPI_MAX,0,IssmComm::GetComm() );
-        ISSM_MPI_Bcast(&node_minvel,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());   
+        ISSM_MPI_Bcast(&node_minvel,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
         minvel=node_minvel;
 
@@ -1526 +1526 @@
         /*Figure out minimum across the cluster: */
         ISSM_MPI_Reduce(&minvx,&node_minvx,1,ISSM_MPI_DOUBLE,ISSM_MPI_MAX,0,IssmComm::GetComm() );
-        ISSM_MPI_Bcast(&node_minvx,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());   
+        ISSM_MPI_Bcast(&node_minvx,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
         minvx=node_minvx;
 
@@ -1550 +1550 @@
         /*Figure out minimum across the cluster: */
         ISSM_MPI_Reduce(&minvy,&node_minvy,1,ISSM_MPI_DOUBLE,ISSM_MPI_MAX,0,IssmComm::GetComm() );
-        ISSM_MPI_Bcast(&node_minvy,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());   
+        ISSM_MPI_Bcast(&node_minvy,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
         minvy=node_minvy;
 
@@ -1574 +1574 @@
         /*Figure out minimum across the cluster: */
         ISSM_MPI_Reduce(&minvz,&node_minvz,1,ISSM_MPI_DOUBLE,ISSM_MPI_MAX,0,IssmComm::GetComm() );
-        ISSM_MPI_Bcast(&node_minvz,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());   
+        ISSM_MPI_Bcast(&node_minvz,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
         minvz=node_minvz;
 
@@ -1619 +1619 @@
                         omega_input->GetInputDerivativeValue(&dp[0],xyz_list,gauss);
 
-                        /*Tikhonov regularization: J = 1/2 ((dp/dx)^2 + (dp/dy)^2) */ 
+                        /*Tikhonov regularization: J = 1/2 ((dp/dx)^2 + (dp/dy)^2) */
                         //J+=weight*1/2*(dp[0]*dp[0]+dp[1]*dp[1])*Jdet*gauss->weight;
                         J+=weight*1/2*pow(dp[0]*dp[0]+dp[1]*dp[1],2)*Jdet*gauss->weight;
@@ -1678 +1678 @@
                         omega0_input->GetInputValue(&p0,gauss);
 
-                        /*Tikhonov regularization: J = 1/2 ((dp/dx)^2 + (dp/dy)^2) */ 
+                        /*Tikhonov regularization: J = 1/2 ((dp/dx)^2 + (dp/dy)^2) */
                         //J+=weight*1/2*(dp[0]*dp[0]+dp[1]*dp[1])*Jdet*gauss->weight;
                         J+=weight*1/2*pow(p - p0,2)*Jdet*gauss->weight;
@@ -1902 +1902 @@
                                                 ISSM_MPI_Bcast(&nodesperelement,1,ISSM_MPI_INT,0,IssmComm::GetComm());
                                                 ISSM_MPI_Bcast(&array_size,1,ISSM_MPI_INT,0,IssmComm::GetComm());
-                                                
+
                                                 if(results_on_nodes){
 
@@ -1955 +1955 @@
                                                                 IssmDouble* values    = xNewZeroInit<IssmDouble>(nlines*ncols);
                                                                 IssmDouble* allvalues = xNew<IssmDouble>(nlines*ncols);
-                                                                
+
                                                                 /*Fill-in matrix*/
                                                                 for(int j=0;j<elements->Size();j++){
@@ -1964 +1964 @@
                                                                 ISSM_MPI_Allreduce((void*)values,(void*)allvalues,ncols*nlines,ISSM_MPI_PDOUBLE,ISSM_MPI_SUM,IssmComm::GetComm());
                                                                 xDelete<IssmDouble>(values);
-                                                                
+
                                                                 if(save_results)results->AddResult(new GenericExternalResult<IssmDouble*>(results->Size()+1,output_enum,allvalues,nlines,ncols,step,time));
                                                                 xDelete<IssmDouble>(allvalues);
@@ -2012 +2012 @@
         int domaintype;
         this->parameters->FindParam(&domaintype,DomainTypeEnum);
-        
+
         /*1: go throug all elements of this partition and figure out how many
          * segments we have (corresopnding to levelset = 0)*/
@@ -2132 +2132 @@
                 case VelEnum:                            this->ElementResponsex(responses,VelEnum); break;
                 case FrictionCoefficientEnum:            NodalValuex(responses, FrictionCoefficientEnum,elements,nodes, vertices, loads, materials, parameters); break;
-                default: 
+                default:
                         if(response_descriptor_enum>=Outputdefinition1Enum && response_descriptor_enum <=Outputdefinition100Enum){
                                 IssmDouble double_result = OutputDefinitionsResponsex(this,response_descriptor_enum);
                                 *responses=double_result;
                         }
-                        else _error_("response descriptor \"" << EnumToStringx(response_descriptor_enum) << "\" not supported yet!"); 
-                        break; 
+                        else _error_("response descriptor \"" << EnumToStringx(response_descriptor_enum) << "\" not supported yet!");
+                        break;
         }
 
@@ -2269 +2269 @@
                         thickness_input->GetInputDerivativeValue(&dp[0],xyz_list,gauss);
 
-                        /*Tikhonov regularization: J = 1/2 ((dp/dx)^2 + (dp/dy)^2) */ 
+                        /*Tikhonov regularization: J = 1/2 ((dp/dx)^2 + (dp/dy)^2) */
                         J+=weight*1/2*(dp[0]*dp[0]+dp[1]*dp[1])*Jdet*gauss->weight;
                 }
@@ -2460 +2460 @@
         analysis->UpdateConstraints(this);
         delete analysis;
-        
+
         /*Second, constraints might be time dependent: */
-        SpcNodesx(nodes,constraints,parameters,analysis_type); 
+        SpcNodesx(nodes,constraints,parameters,analysis_type);
 
         /*Now, update degrees of freedoms: */
@@ -2473 +2473 @@
         IssmDouble         *surface = NULL;
         IssmDouble         *bed     = NULL;
-                        
+
         if(VerboseSolution()) _printf0_("   updating vertices positions\n");
 
@@ -2512 +2512 @@
 
 /*AMR*/
-#if !defined(_HAVE_ADOLC_) 
+#if !defined(_HAVE_ADOLC_)
 void FemModel::ReMesh(void){/*{{{*/
 
@@ -2522 +2522 @@
         int newnumberofvertices = -1;
         int newnumberofelements = -1;
-        bool* my_elements                       = NULL; 
+        bool* my_elements                       = NULL;
         int* my_vertices                        = NULL;
         int elementswidth       = this->GetElementsWidth();//just tria elements in this version
@@ -2528 +2528 @@
         bool isgroundingline;
 
-        /*Branch to specific amr depending on requested method*/        
+        /*Branch to specific amr depending on requested method*/
         this->parameters->FindParam(&amrtype,AmrTypeEnum);
         switch(amrtype){
@@ -2541 +2541 @@
                 default: _error_("not implemented yet");
         }
-        
+
         /*Partitioning the new mesh. Maybe ElementsAndVerticesPartitioning.cpp could be modified to set this without iomodel.*/
         this->ElementsAndVerticesPartitioning(newnumberofvertices,newnumberofelements,elementswidth,newelementslist,&my_elements,&my_vertices);
@@ -2572 +2572 @@
 
                 /*As the domain is 2D, it is not necessary to create nodes for this analysis*/
-                if(analysis_enum==StressbalanceVerticalAnalysisEnum) continue;    
+                if(analysis_enum==StressbalanceVerticalAnalysisEnum) continue;
 
                 this->CreateNodes(newnumberofvertices,my_vertices,nodecounter,analysis_enum,new_nodes);
@@ -2602 +2602 @@
                 //SetCurrentConfiguration(analysis_type);
 
-                this->analysis_counter=i;       
+                this->analysis_counter=i;
                 /*Now, plug analysis_counter and analysis_type inside the parameters: */
                 this->parameters->SetParam(this->analysis_counter,AnalysisCounterEnum);
@@ -2619 +2619 @@
 
                 ConfigureObjectsx(new_elements,this->loads,new_nodes,new_vertices,new_materials,this->parameters);
-                if(i==0){ 
+                if(i==0){
                         VerticesDofx(new_vertices,this->parameters); //only call once, we only have one set of vertices
                 }
@@ -2663 +2663 @@
         /*Insert bedrock from mismip+ setup*/
         /*This was used to Misomip project/simulations*/
-        
+
         if(VerboseSolution())_printf0_("        call Mismip bedrock adjust module\n");
 
@@ -2680 +2680 @@
                         by              = 500./(1.+exp((-2./4000.)*(y-80000./2.-24000.)))+500./(1.+exp((2./4000.)*(y-80000./2.+24000.)));
                         r[i]    = max(bx+by,-720.);
-                }       
+                }
                 /*insert new bedrock*/
                 element->AddInput(BedEnum,&r[0],P1Enum);
@@ -2693 +2693 @@
 
         if(VerboseSolution())_printf0_("        call adjust base and thickness module\n");
-        
+
         int     numvertices = this->GetElementsWidth();
    IssmDouble rho_water,rho_ice,density,base_float;
@@ -2705 +2705 @@
         for(int el=0;el<this->elements->Size();el++){
                 Element* element=xDynamicCast<Element*>(this->elements->GetObjectByOffset(el));
-        
+
                 element->GetInputListOnVertices(&s[0],SurfaceEnum);
                 element->GetInputListOnVertices(&r[0],BedEnum);
@@ -2717 +2717 @@
                         base_float = rho_ice*s[i]/(rho_ice-rho_water);
                         if(r[i]>base_float){
-                                b[i] = r[i];                    
-                        } 
+                                b[i] = r[i];
+                        }
                         else {
-                                b[i] = base_float;      
-                        } 
+                                b[i] = base_float;
+                        }
 
                         if(fabs(sl[i])>0) _error_("Sea level value "<<sl[i]<<" not supported!");
                         /*update thickness and mask grounded ice level set*/
                         h[i]      = s[i]-b[i];
-                        phi[i]  = h[i]+r[i]/density;    
+                        phi[i]  = h[i]+r[i]/density;
                 }
 
@@ -2734 +2734 @@
                 element->AddInput(BaseEnum,&b[0],P1Enum);
         }
-        
+
    /*Delete*/
    xDelete<IssmDouble>(phi);
@@ -2762 +2762 @@
         Vector<IssmDouble>* input_interpolations        = NULL;
         IssmDouble* input_interpolations_serial = NULL;
-   int* pos                                                                                             = NULL; 
+   int* pos                                                                                             = NULL;
         IssmDouble value                                                                        = 0;
 
@@ -2782 +2782 @@
                         P0input_interp = xNew<int>(numP0inputs);
                         P1input_enums  = xNew<int>(numP1inputs);
-                        P1input_interp = xNew<int>(numP1inputs);        
+                        P1input_interp = xNew<int>(numP1inputs);
                 }
                 numP0inputs = 0;
@@ -2814 +2814 @@
                 }
         }
-        
-        /*Get P0 and P1 inputs over the elements*/      
+
+        /*Get P0 and P1 inputs over the elements*/
         pos             = xNew<int>(elementswidth);
         vP0inputs= new Vector<IssmDouble>(numberofelements*numP0inputs);
@@ -2821 +2821 @@
         for(int i=0;i<this->elements->Size();i++){
                 Element* element=xDynamicCast<Element*>(this->elements->GetObjectByOffset(i));
-                
+
                 /*Get P0 inputs*/
                 for(int j=0;j<numP0inputs;j++){
-                        TriaInput* input=xDynamicCast<TriaInput*>(element->GetInput(P0input_enums[j]));         
+                        TriaInput* input=xDynamicCast<TriaInput*>(element->GetInput(P0input_enums[j]));
                         input->GetInputAverage(&value);
                         pos[0]=element->Sid()*numP0inputs+j;
-                        /*Insert input in the vector*/  
+                        /*Insert input in the vector*/
                         vP0inputs->SetValues(1,pos,&value,INS_VAL);
                 }
-                
+
                 /*Get P1 inputs*/
                 for(int j=0;j<numP1inputs;j++){
@@ -2837 +2837 @@
                         pos[1]=element->vertices[1]->Sid()*numP1inputs+j;
                         pos[2]=element->vertices[2]->Sid()*numP1inputs+j;
-                        /*Insert input in the vector*/  
-                        vP1inputs->SetValues(elementswidth,pos,input->values,INS_VAL);  
+                        /*Insert input in the vector*/
+                        vP1inputs->SetValues(elementswidth,pos,input->values,INS_VAL);
                 }
         }
@@ -2847 +2847 @@
         P0inputs=vP0inputs->ToMPISerial();
         P1inputs=vP1inputs->ToMPISerial();
-        
+
         /*Assign pointers*/
-        *pnumP0inputs           = numP0inputs; 
-        *pP0inputs                      = P0inputs; 
+        *pnumP0inputs           = numP0inputs;
+        *pP0inputs                      = P0inputs;
         *pP0input_enums = P0input_enums;
         *pP0input_interp        = P0input_interp;
-        *pnumP1inputs           = numP1inputs; 
-        *pP1inputs                      = P1inputs; 
+        *pnumP1inputs           = numP1inputs;
+        *pP1inputs                      = P1inputs;
         *pP1input_enums = P1input_enums;
-        *pP1input_interp        = P1input_interp;       
+        *pP1input_interp        = P1input_interp;
 
         /*Cleanup*/
@@ -2867 +2867 @@
 /*}}}*/
 void FemModel::InterpolateInputs(Vertices* newfemmodel_vertices,Elements* newfemmodel_elements){/*{{{*/
-        
+
         int numberofelements                    = this->elements->NumberOfElements();   //global, entire old mesh
         int newnumberofelements         = newfemmodel_elements->Size();                 //just on the new partition
@@ -2883 +2883 @@
         int* P1input_enums                      = NULL;
         int* P1input_interp                     = NULL;
-        IssmDouble* values                      = NULL; 
+        IssmDouble* values                      = NULL;
    IssmDouble* vector           = NULL;
         IssmDouble* x                                   = NULL;//global, entire old mesh
@@ -2895 +2895 @@
         IssmDouble* newyc                               = NULL;//just on the new partition
         int* newelementslist                    = NULL;//just on the new partition
-        int* sidtoindex                         = NULL;//global vertices sid to partition index 
+        int* sidtoindex                         = NULL;//global vertices sid to partition index
 
         /*Get old P0 and P1  inputs (entire mesh)*/
@@ -2928 +2928 @@
 
         /*Insert P0 and P1 inputs into the new elements (just on the new partition)*/
-        values=xNew<IssmDouble>(elementswidth); 
+        values=xNew<IssmDouble>(elementswidth);
         for(int i=0;i<newfemmodel_elements->Size();i++){//just on the new partition
                 Element* element=xDynamicCast<Element*>(newfemmodel_elements->GetObjectByOffset(i));
                 /*newP0inputs is just on the new partition*/
                 for(int j=0;j<numP0inputs;j++){
-                        switch(P0input_interp[j]){      
+                        switch(P0input_interp[j]){
                                 case P0Enum:
                                 case DoubleInputEnum:
                                         element->AddInput(new DoubleInput(P0input_enums[j],newP0inputs[i*numP0inputs+j]));
                                         break;
-                                case IntInputEnum: 
+                                case IntInputEnum:
                                         element->AddInput(new IntInput(P0input_enums[j],reCast<int>(newP0inputs[i*numP0inputs+j])));
                                         break;
@@ -2956 +2956 @@
                 }
         }
-        
+
         /*Cleanup*/
         xDelete<IssmDouble>(P0inputs);
@@ -2995 +2995 @@
 
         if(!this->elements || !this->vertices || !this->results || !this->parameters) return;
-         
+
         parameters->FindParam(&step,StepEnum);
         parameters->FindParam(&time,TimeEnum);
@@ -3013 +3013 @@
         this->results->AddResult(new GenericExternalResult<IssmDouble*>(this->results->Size()+1,MeshYEnum,
                                         y,numberofvertices,1,step,time));
-        
+
         /*Cleanup*/
         xDelete<IssmDouble>(x);
@@ -3074 +3074 @@
    /*Assemble*/
         vmasklevelset->Assemble();
-        
+
         /*Serialize and set output*/
         (*pmasklevelset)=vmasklevelset->ToMPISerial();
@@ -3088 +3088 @@
 
         /*newelementslist is in Matlab indexing*/
-        
+
         /*Creating connectivity table*/
         int* connectivity=NULL;
@@ -3099 +3099 @@
                         connectivity[vertexid-1]+=1;//Matlab to C indexing
                 }
-        }       
+        }
 
         /*Create vertex and insert in vertices*/
         for(int i=0;i<newnumberofvertices;i++){
                 if(my_vertices[i]){
-                        Vertex *newvertex=new Vertex(); 
+                        Vertex *newvertex=new Vertex();
                         newvertex->id=i+1;
                         newvertex->sid=i;
@@ -3115 +3115 @@
                         newvertex->connectivity=connectivity[i];
                         newvertex->clone=false;//itapopo check this
-                        vertices->AddObject(newvertex); 
-                } 
+                        vertices->AddObject(newvertex);
+                }
         }
 
@@ -3143 +3143 @@
                         }
                         else newtria->element_type_list=NULL;
-                        
+
                         /*Element hook*/
                         int matpar_id=newnumberofelements+1; //retrieve material parameter id (last pointer in femodel->materials)
@@ -3149 +3149 @@
                         /*retrieve vertices ids*/
                         int* vertex_ids=xNew<int>(elementswidth);
-                        for(int j=0;j<elementswidth;j++)        vertex_ids[j]=reCast<int>(newelementslist[elementswidth*i+j]);//this Hook wants Matlab indexing 
+                        for(int j=0;j<elementswidth;j++)        vertex_ids[j]=reCast<int>(newelementslist[elementswidth*i+j]);//this Hook wants Matlab indexing
                         /*Setting the hooks*/
                         newtria->numanalyses =this->nummodels;
@@ -3161 +3161 @@
                         /*Clean up*/
                         xDelete<int>(vertex_ids);
-                        elements->AddObject(newtria);   
-                } 
+                        elements->AddObject(newtria);
+                }
         }
 
@@ -3172 +3172 @@
         for(int i=0;i<newnumberofelements;i++){
                 if(my_elements[i]){
-                        materials->AddObject(new Matice(i+1,i,MaticeEnum));     
-                } 
-        }
-        
+                        materials->AddObject(new Matice(i+1,i,MaticeEnum));
+                }
+        }
+
         /*Add new constant material property to materials, at the end: */
         Matpar *newmatpar=static_cast<Matpar*>(this->materials->GetObjectByOffset(this->materials->Size()-1)->copy());
         newmatpar->SetMid(newnumberofelements+1);
-        materials->AddObject(newmatpar);//put it at the end of the materials        
+        materials->AddObject(newmatpar);//put it at the end of the materials
 }
 /*}}}*/
@@ -3186 +3186 @@
         int lid=0;
         for(int j=0;j<newnumberofvertices;j++){
-                if(my_vertices[j]){                             
-                        
-                        Node* newnode=new Node();       
-                        
+                if(my_vertices[j]){
+
+                        Node* newnode=new Node();
+
                         /*id: */
                         newnode->id=nodecounter+j+1;
@@ -3195 +3195 @@
                         newnode->lid=lid++;
                         newnode->analysis_enum=analysis_enum;
-                        
+
                         /*Initialize coord_system: Identity matrix by default*/
                         for(int k=0;k<3;k++) for(int l=0;l<3;l++) newnode->coord_system[k][l]=0.0;
                         for(int k=0;k<3;k++) newnode->coord_system[k][k]=1.0;
-                        
+
                         /*indexing:*/
                         newnode->indexingupdate=true;
-                        
+
                         Analysis* analysis=EnumToAnalysis(analysis_enum);
                         int *doftypes=NULL;
@@ -3216 +3216 @@
                         /*Stressbalance Horiz*/
                         if(analysis_enum==StressbalanceAnalysisEnum){
-                                // itapopo this code is rarely used. 
+                                // itapopo this code is rarely used.
                                 /*Coordinate system provided, convert to coord_system matrix*/
                                 //XZvectorsToCoordinateSystem(&this->coord_system[0][0],&iomodel->Data(StressbalanceReferentialEnum)[j*6]);
@@ -3231 +3231 @@
         if(!femmodel_vertices) _error_("GetMesh: vertices are NULL.");
         if(!femmodel_elements) _error_("GetMesh: elements are NULL.");
-        
+
         int numberofvertices = femmodel_vertices->NumberOfVertices();
         int numberofelements = femmodel_elements->NumberOfElements();
@@ -3237 +3237 @@
         IssmDouble* x                   = NULL;
         IssmDouble* y                   = NULL;
-        IssmDouble* z                   = NULL; 
+        IssmDouble* z                   = NULL;
         int* elementslist       = NULL;
         int* elem_vertices      = NULL;
@@ -3246 +3246 @@
         /*Get vertices coordinates*/
         VertexCoordinatesx(&x,&y,&z,femmodel_vertices,false) ;
-        
+
         /*Get element vertices*/
         elem_vertices                           = xNew<int>(elementswidth);
@@ -3261 +3261 @@
         vid3->SetValue(element->sid,elem_vertices[2],INS_VAL);
    }
-                
+
         /*Assemble*/
    vid1->Assemble();
@@ -3271 +3271 @@
    id2 = vid2->ToMPISerial();
         id3 = vid3->ToMPISerial();
-        
+
         /*Construct elements list*/
         elementslist=xNew<int>(numberofelements*elementswidth);
@@ -3280 +3280 @@
                 elementslist[elementswidth*i+2] = reCast<int>(id3[i])+1; //InterpMesh wants Matlab indexing
         }
-        
+
         /*Assign pointers*/
         *px                             = x;
@@ -3301 +3301 @@
         if(!femmodel_vertices) _error_("GetMesh: vertices are NULL.");
         if(!femmodel_elements) _error_("GetMesh: elements are NULL.");
-        
+
         int numberofvertices                    = femmodel_vertices->Size();    //number of vertices of this partition
         int numbertotalofvertices       = femmodel_vertices->NumberOfVertices();        //number total of vertices (entire mesh)
@@ -3308 +3308 @@
         IssmDouble* x                                   = NULL;
         IssmDouble* y                                   = NULL;
-        IssmDouble* z                                   = NULL; 
+        IssmDouble* z                                   = NULL;
         int* elementslist                               = NULL;
         int* sidtoindex                         = NULL;
         int* elem_vertices                      = NULL;
-        
+
         /*Get vertices coordinates of this partition*/
         sidtoindex      = xNewZeroInit<int>(numbertotalofvertices);//entire mesh, all vertices
@@ -3318 +3318 @@
         y                               = xNew<IssmDouble>(numberofvertices);//just this partitio;
         z                               = xNew<IssmDouble>(numberofvertices);//just this partitio;
-        
+
         /*Go through in this partition (vertices)*/
         for(int i=0;i<numberofvertices;i++){//just this partition
-                Vertex* vertex=(Vertex*)femmodel_vertices->GetObjectByOffset(i);        
+                Vertex* vertex=(Vertex*)femmodel_vertices->GetObjectByOffset(i);
                 /*Attention: no spherical coordinates*/
                 x[i]=vertex->GetX();
@@ -3334 +3334 @@
         elementslist = xNew<int>(numberofelements*elementswidth);
         if(numberofelements*elementswidth<0) _error_("numberofelements negative.");
-        
+
         for(int i=0;i<numberofelements;i++){//just this partition
         Element* element=xDynamicCast<Element*>(femmodel_elements->GetObjectByOffset(i));
@@ -3341 +3341 @@
                 elementslist[elementswidth*i+1] = sidtoindex[elem_vertices[1]]+1; //InterpMesh wants Matlab indexing
                 elementslist[elementswidth*i+2] = sidtoindex[elem_vertices[2]]+1; //InterpMesh wants Matlab indexing
-        }       
-                
+        }
+
         /*Assign pointers*/
         *px                             = x;
@@ -3348 +3348 @@
         *pz                             = z;
         *pelementslist = elementslist; //Matlab indexing. InterMesh uses this type.
-        *psidtoindex    = sidtoindex;  //it is ncessary to insert inputs 
+        *psidtoindex    = sidtoindex;  //it is ncessary to insert inputs
 
         /*Cleanup*/
@@ -3359 +3359 @@
         /*OTHERS CONSTRAINTS MUST BE IMPLEMENTED*/
         if(analysis_enum!=StressbalanceAnalysisEnum) return;
-        
+
         int numberofnodes_analysistype= this->nodes->NumberOfNodes(analysis_enum);
-        int dofpernode                                          = 2;                                                                                                            //vx and vy 
+        int dofpernode                                          = 2;                                                                                                            //vx and vy
         int numberofcols                                        = dofpernode*2;                                                                         //to keep dofs and flags in the vspc vector
         int numberofvertices                            = this->vertices->NumberOfVertices();                   //global, entire old mesh
@@ -3385 +3385 @@
         newy=xNew<IssmDouble>(newnumberofvertices);//just the new partition
         for(int i=0;i<newnumberofvertices;i++){//just the new partition
-                Vertex* vertex=(Vertex*)newfemmodel_vertices->GetObjectByOffset(i);     
+                Vertex* vertex=(Vertex*)newfemmodel_vertices->GetObjectByOffset(i);
                 /*Attention: no spherical coordinates*/
                 newx[i]=vertex->GetX();
@@ -3393 +3393 @@
         /*Get spcvx and spcvy of old mesh*/
         for(int i=0;i<this->constraints->Size();i++){
-                
+
                 Constraint* constraint=(Constraint*)constraints->GetObjectByOffset(i);
                 if(!constraint->InAnalysis(analysis_enum)) _error_("AMR create constraints for "<<EnumToStringx(analysis_enum)<<" not supported yet!\n");
@@ -3400 +3400 @@
                 int dof                                 = spcstatic->GetDof();
                 int node                                        = spcstatic->GetNodeId();
-                IssmDouble spcvalue     = spcstatic->GetValue(); 
+                IssmDouble spcvalue     = spcstatic->GetValue();
                 int nodeindex                   = node-1;
-                
+
                 /*vx and vx flag insertion*/
                 if(dof==0) {//vx
                         vspc->SetValue(nodeindex*numberofcols,spcvalue,INS_VAL);    //vx
                         vspc->SetValue(nodeindex*numberofcols+dofpernode,1,INS_VAL);//vxflag
-                } 
+                }
                 /*vy and vy flag insertion*/
                 if(dof==1){//vy
@@ -3423 +3423 @@
                                                                 spc,numberofvertices,numberofcols,
                                                                 newx,newy,newnumberofvertices,NULL);
-        
+
         /*Now, insert the interpolated constraints in the data set (constraints)*/
         count=0;
@@ -3440 +3440 @@
                 /*spcvy*/
                 if(!xIsNan<IssmDouble>(newspc[i*numberofcols+1]) && newspc[i*numberofcols+dofpernode+1]>(1-eps) ){
-                        newfemmodel_constraints->AddObject(new SpcStatic(constraintcounter+count+1,nodecounter+vertex->Sid()+1,1,newspc[i*numberofcols+1],analysis_enum)); 
+                        newfemmodel_constraints->AddObject(new SpcStatic(constraintcounter+count+1,nodecounter+vertex->Sid()+1,1,newspc[i*numberofcols+1],analysis_enum));
                         //add count'th spc, on node i+1, setting dof 1 to vx.
                         count++;
@@ -3497 +3497 @@
         bool *my_elements = NULL;
         int *my_vertices  = NULL;
-        
-        _assert_(newnumberofvertices>0); 
-        _assert_(newnumberofelements>0); 
+
+        _assert_(newnumberofvertices>0);
+        _assert_(newnumberofelements>0);
         epart=xNew<int>(newnumberofelements);
         npart=xNew<int>(newnumberofvertices);
    index=xNew<int>(elementswidth*newnumberofelements);
-   
+
         for (int i=0;i<newnumberofelements;i++){
         for (int j=0;j<elementswidth;j++){
@@ -3523 +3523 @@
                 for (int i=0;i<newnumberofvertices;i++) npart[i]=0;
         }
-        else _error_("At least one processor is required");         
+        else _error_("At least one processor is required");
 
         my_vertices=xNew<int>(newnumberofvertices);
@@ -3533 +3533 @@
         for(int i=0;i<newnumberofelements;i++){
                 /*!All elements have been partitioned above, only deal with elements for this cpu: */
-                if(my_rank==epart[i]){ 
+                if(my_rank==epart[i]){
                         my_elements[i]=true;
-                        /*Now that we are here, we can also start building the list of vertices belonging to this cpu partition: we use 
-                         *the  element index to do this. For each element n, we know index[n][0:2] holds the indices (matlab indexing) 
-                         into the vertices coordinates. If we start plugging 1 into my_vertices for each index[n][i] (i=0:2), then my_vertices 
+                        /*Now that we are here, we can also start building the list of vertices belonging to this cpu partition: we use
+                         *the  element index to do this. For each element n, we know index[n][0:2] holds the indices (matlab indexing)
+                         into the vertices coordinates. If we start plugging 1 into my_vertices for each index[n][i] (i=0:2), then my_vertices
                          will hold which vertices belong to this partition*/
                         for(int j=0;j<elementswidth;j++){
-                                _assert_(newelementslist[elementswidth*i+j]-1<newnumberofvertices);//newelementslist is in Matlab indexing 
+                                _assert_(newelementslist[elementswidth*i+j]-1<newnumberofvertices);//newelementslist is in Matlab indexing
                                 my_vertices[newelementslist[elementswidth*i+j]-1]=1;//newelementslist is in Matlab indexing
                         }
@@ -3552 +3552 @@
         /*Free ressources:*/
         xDelete<int>(epart);
-        xDelete<int>(npart);        
+        xDelete<int>(npart);
         xDelete<int>(index);
 }
 /*}}}*/
 void FemModel::SmoothedDeviatoricStressTensor(IssmDouble** ptauxx,IssmDouble** ptauyy,IssmDouble** ptauxy){/*{{{*/
-        
+
         int elementswidth                                                       = this->GetElementsWidth();//just 2D mesh, tria elements
    int numberofvertices                                         = this->vertices->NumberOfVertices();
    IssmDouble weight                                            = 0.;
-        IssmDouble*     tauxx                                                   = NULL; 
-        IssmDouble*     tauyy                                                   = NULL; 
-        IssmDouble*     tauxy                                                   = NULL; 
+        IssmDouble*     tauxx                                                   = NULL;
+        IssmDouble*     tauyy                                                   = NULL;
+        IssmDouble*     tauxy                                                   = NULL;
    IssmDouble* totalweight                              = NULL;
         IssmDouble* deviatoricstressxx          = xNew<IssmDouble>(elementswidth);
@@ -3573 +3573 @@
    Vector<IssmDouble>* vectauxy                 = new Vector<IssmDouble>(numberofvertices);
    Vector<IssmDouble>* vectotalweight   = new Vector<IssmDouble>(numberofvertices);
-        
+
         /*Update the Deviatoric Stress tensor over the elements*/
         this->DeviatoricStressx();
-        
+
    /*Calculate the Smoothed Deviatoric Stress tensor*/
         for(int i=0;i<this->elements->Size();i++){
@@ -3621 +3621 @@
         /*Divide for the total weight*/
         for(int i=0;i<numberofvertices;i++){
-                _assert_(totalweight[i]>0);     
+                _assert_(totalweight[i]>0);
                 tauxx[i] = tauxx[i]/totalweight[i];
                 tauyy[i] = tauyy[i]/totalweight[i];
@@ -3646 +3646 @@
 void FemModel::ZZErrorEstimator(IssmDouble** pelementerror){/*{{{*/
 
-        /*Compute the Zienkiewicz and Zhu (ZZ) error estimator for the deviatoric stress tensor. 
+        /*Compute the Zienkiewicz and Zhu (ZZ) error estimator for the deviatoric stress tensor.
          * Ref.: Zienkiewicz and Zhu, A Simple Error Estimator and Adaptive Procedure for Practical Engineering Analysis, Int. J. Numer. Meth. Eng, 1987*/
 
@@ -3666 +3666 @@
         /*Get smoothed deviatoric stress tensor*/
         this->SmoothedDeviatoricStressTensor(&smoothedtauxx,&smoothedtauyy,&smoothedtauxy);
-        
+
         /*Integrate the error over elements*/
    for(int i=0;i<this->elements->Size();i++){
@@ -3674 +3674 @@
       element->GetInputListOnVertices(tauxy,DeviatoricStressxyEnum);
       element->GetVerticesSidList(elem_vertices);
-                
+
                 /*Integrate*/
                 element->GetVerticesCoordinates(&xyz_list);
@@ -3689 +3689 @@
                                 ftxy+=(tauxy[n]-smoothedtauxy[elem_vertices[n]])*basis[n];
                         }
-                        error+=Jdet*gauss->weight*( pow(ftxx,2)+pow(ftyy,2)+pow(ftxy,2) ); //e^2 
-                }
-                /*Set the error in the global vector*/  
+                        error+=Jdet*gauss->weight*( pow(ftxx,2)+pow(ftyy,2)+pow(ftxy,2) ); //e^2
+                }
+                /*Set the error in the global vector*/
       sid=element->Sid();
                 error = sqrt(error);//sqrt(e^2)
@@ -3705 +3705 @@
    /*Serialize and set output*/
    (*pelementerror)=velementerror->ToMPISerial();
-        
+
         /*Cleanup*/
         xDelete<IssmDouble>(smoothedtauxx);
@@ -3749 +3749 @@
       Tria* triaelement = xDynamicCast<Tria*>(element);
       weight            = triaelement->GetArea();//the tria area is a choice for the weight
-      
+
                 /*dH/dx*/
       vecdHdx->SetValue(elem_vertices[0],weight*GradH[0],ADD_VAL);
@@ -3817 +3817 @@
    /*Get smoothed deviatoric stress tensor*/
    this->SmoothedGradThickness(&smoothed_dHdx,&smoothed_dHdy);
-   
+
         /*Integrate the error over elements*/
    for(int i=0;i<this->elements->Size();i++){
@@ -3903 +3903 @@
         IssmDouble* x                                   = NULL;
         IssmDouble* y                                   = NULL;
-        IssmDouble* z                                   = NULL; 
+        IssmDouble* z                                   = NULL;
         IssmDouble* xyz_list                    = NULL;
         IssmDouble x1,y1,x2,y2,x3,y3;
@@ -3912 +3912 @@
       //element->GetVerticesSidList(elem_vertices);
       int sid = element->Sid();
-                element->GetVerticesCoordinates(&xyz_list); 
+                element->GetVerticesCoordinates(&xyz_list);
                 x1 = xyz_list[3*0+0];y1 = xyz_list[3*0+1];
                 x2 = xyz_list[3*1+0];y2 = xyz_list[3*1+1];
@@ -3945 +3945 @@
                 _error_("level set type not implemented yet!");
         }
-        
+
         /*Outputs*/
         IssmDouble* zerolevelset_points                 = NULL;
         int npoints                                                                             = 0;
-        
+
         /*Intermediaries*/
         int elementswidth                       = this->GetElementsWidth();
@@ -3962 +3962 @@
         int count,sid;
         IssmDouble xc,yc,x1,y1,x2,y2,x3,y3;
-        
+
         /*Use the element center coordinate if level set is zero (grounding line or ice front), otherwise set NAN*/
    for(int i=0;i<this->elements->Size();i++){
@@ -3969 +3969 @@
                 element->GetVerticesSidList(elem_vertices);
                 sid= element->Sid();
-                element->GetVerticesCoordinates(&xyz_list); 
+                element->GetVerticesCoordinates(&xyz_list);
                 x1 = xyz_list[3*0+0];y1 = xyz_list[3*0+1];
                 x2 = xyz_list[3*1+0];y2 = xyz_list[3*1+1];
                 x3 = xyz_list[3*2+0];y3 = xyz_list[3*2+1];
                 xc      = NAN;
-                yc      = NAN;  
+                yc      = NAN;
         Tria* tria      = xDynamicCast<Tria*>(element);
                 if(tria->IsIceInElement()){/*verify if there is ice in the element*/
-                        if(levelset[0]*levelset[1]<0. || levelset[0]*levelset[2]<0. ||  
+                        if(levelset[0]*levelset[1]<0. || levelset[0]*levelset[2]<0. ||
                                 abs(levelset[0]*levelset[1])<DBL_EPSILON || abs(levelset[0]*levelset[2])<DBL_EPSILON) {
                                 xc=(x1+x2+x3)/3.;
@@ -4007 +4007 @@
                 }
         }
-        
+
         /*Assign outputs*/
         numberofpoints                          = npoints;
@@ -4047 +4047 @@
         responses_pointer=d_responses;
 
-        //watch out, we have more d_numresponses than numresponsedescriptors, because the responses have been expanded if they were scaled. 
+        //watch out, we have more d_numresponses than numresponsedescriptors, because the responses have been expanded if they were scaled.
         //because we don't know the d_responses descriptors (the scaled ones) we can't key off them, so we will key off the responses_descriptors: */
 
@@ -4140 +4140 @@
 
         int         ns,nsmax;
-        
+
         /*Go through elements, and add contribution from each element to the deflection vector wg:*/
         ns = elements->Size();
-        
+
         /*Figure out max of ns: */
         ISSM_MPI_Reduce(&ns,&nsmax,1,ISSM_MPI_INT,ISSM_MPI_MAX,0,IssmComm::GetComm());
@@ -4162 +4162 @@
                 }
         }
-        
+
         /*One last time: */
         pUp->Assemble();
@@ -4181 +4181 @@
 
         int         ns,nsmax;
-        
+
         /*Go through elements, and add contribution from each element to the deflection vector wg:*/
         ns = elements->Size();
-        
-        /*First, figure out the surface area of Earth: */ 
+
+        /*First, figure out the surface area of Earth: */
         for(int i=0;i<ns;i++){
                 Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
@@ -4209 +4209 @@
                 }
         }
-        
+
         /*One last time: */
         pUp->Assemble();
@@ -4237 +4237 @@
         IssmDouble  eartharea_cpu  = 0.;
         int         ns,nsmax;
-        
+
         /*Go through elements, and add contribution from each element to the deflection vector wg:*/
         ns = elements->Size();
@@ -4261 +4261 @@
         for(int i=0;i<nsmax;i++){
                 if(i<ns){
-                
+
                         if(VerboseConvergence())if(i%100==0)_printf0_("\r" << "      convolution progress: " << (double)i/(double)ns*100 << "%  ");
-                
+
                         Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
                         element->SealevelriseEustatic(pSgi,&eustatic_cpu_e,latitude,longitude,radius,oceanarea,eartharea);
@@ -4271 +4271 @@
         }
         if(VerboseConvergence())_printf0_("\n");
-                
+
         /*One last time: */
         pSgi->Assemble();
@@ -4289 +4289 @@
         /*serialized vectors:*/
         IssmDouble* Sg_old=NULL;
-        
+
         IssmDouble  eartharea=0;
         IssmDouble  eartharea_cpu=0;
 
         int         ns,nsmax;
-        
+
         /*Serialize vectors from previous iteration:*/
         Sg_old=pSg_old->ToMPISerial();
@@ -4306 +4306 @@
                 eartharea_cpu += element->GetAreaSpherical();
         }
-        
+
         ISSM_MPI_Reduce (&eartharea_cpu,&eartharea,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() );
         ISSM_MPI_Bcast(&eartharea,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
@@ -4324 +4324 @@
         }
         if(verboseconvolution)if(VerboseConvergence())_printf_("\n");
-        
+
         /*Free ressources:*/
         xDelete<IssmDouble>(Sg_old);
@@ -4336 +4336 @@
         IssmDouble  eartharea_cpu=0;
         IssmDouble      tide_love_h, tide_love_k, fluid_love, moi_e, moi_p, omega, g;
-        IssmDouble      load_love_k2 = -0.30922675; //degree 2 load Love number 
-        IssmDouble      m1, m2, m3; 
-        IssmDouble      lati, longi, radi, value; 
+        IssmDouble      load_love_k2 = -0.30922675; //degree 2 load Love number
+        IssmDouble      m1, m2, m3;
+        IssmDouble      lati, longi, radi, value;
 
         /*Serialize vectors from previous iteration:*/
@@ -4351 +4351 @@
         ISSM_MPI_Bcast(&eartharea,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
 
-        IssmDouble moi_list[3]={0,0,0}; 
-        IssmDouble moi_list_cpu[3]={0,0,0}; 
+        IssmDouble moi_list[3]={0,0,0};
+        IssmDouble moi_list_cpu[3]={0,0,0};
         for(int i=0;i<elements->Size();i++){
                 Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
                 element->SealevelriseMomentOfInertia(&moi_list[0],Sg_old,eartharea);
-                moi_list_cpu[0] += moi_list[0]; 
-                moi_list_cpu[1] += moi_list[1]; 
-                moi_list_cpu[2] += moi_list[2]; 
+                moi_list_cpu[0] += moi_list[0];
+                moi_list_cpu[1] += moi_list[1];
+                moi_list_cpu[2] += moi_list[2];
         }
         ISSM_MPI_Reduce (&moi_list_cpu[0],&moi_list[0],1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() );
         ISSM_MPI_Bcast(&moi_list[0],1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
-        //      
+        //
         ISSM_MPI_Reduce (&moi_list_cpu[1],&moi_list[1],1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() );
         ISSM_MPI_Bcast(&moi_list[1],1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
-        //      
+        //
         ISSM_MPI_Reduce (&moi_list_cpu[2],&moi_list[2],1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() );
         ISSM_MPI_Bcast(&moi_list[2],1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
-        
+
         /*pull out some useful parameters: */
         parameters->FindParam(&tide_love_h,SealevelriseTidalLoveHEnum);
@@ -4378 +4378 @@
 
         /*compute perturbation terms for angular velocity vector: */
-        m1 = 1/(1-tide_love_k/fluid_love) * (1+load_love_k2)/(moi_p-moi_e) * moi_list[0]; 
-        m2 = 1/(1-tide_love_k/fluid_love) * (1+load_love_k2)/(moi_p-moi_e) * moi_list[1]; 
-        m3 = -(1+load_love_k2)/moi_p * moi_list[2];     // term associated with fluid number (3-order-of-magnitude smaller) is negelected  
+        m1 = 1/(1-tide_love_k/fluid_love) * (1+load_love_k2)/(moi_p-moi_e) * moi_list[0];
+        m2 = 1/(1-tide_love_k/fluid_love) * (1+load_love_k2)/(moi_p-moi_e) * moi_list[1];
+        m3 = -(1+load_love_k2)/moi_p * moi_list[2];     // term associated with fluid number (3-order-of-magnitude smaller) is negelected
 
         /* Green's function (1+k_2-h_2/g): checked against Glenn Milne's thesis Chapter 3 (eqs: 3.3-4, 3.10-11)
-         * Perturbation terms for angular velocity vector (m1, m2, m3): checked against Mitrovica (2005 Appendix) & Jensen et al (2013 Appendix A3) 
-         * Sea level rotational feedback: checked against GMD eqs 8-9 (only first order terms, i.e., degree 2 order 0 & 1 considered) 
-         * all DONE in Geographic coordinates: theta \in [-90,90], lambda \in [-180 180] 
+         * Perturbation terms for angular velocity vector (m1, m2, m3): checked against Mitrovica (2005 Appendix) & Jensen et al (2013 Appendix A3)
+         * Sea level rotational feedback: checked against GMD eqs 8-9 (only first order terms, i.e., degree 2 order 0 & 1 considered)
+         * all DONE in Geographic coordinates: theta \in [-90,90], lambda \in [-180 180]
          */
         for(int i=0;i<vertices->Size();i++){
@@ -4395 +4395 @@
                 lati=latitude[sid]/180*PI;      longi=longitude[sid]/180*PI; radi=radius[sid];
 
-                /*only first order terms are considered now: */ 
+                /*only first order terms are considered now: */
                 value=((1.0+tide_love_k-tide_love_h)/9.81)*pow(omega*radi,2.0)*
-                                                (-m3/6.0 + 0.5*m3*cos(2.0*lati) - 0.5*sin(2.*lati)*(m1*cos(longi)+m2*sin(longi))); 
-        
+                                                (-m3/6.0 + 0.5*m3*cos(2.0*lati) - 0.5*sin(2.*lati)*(m1*cos(longi)+m2*sin(longi)));
+
                 pSgo_rot->SetValue(sid,value,INS_VAL); //INS_VAL ensures that you don't add several times
         }
@@ -4404 +4404 @@
         /*Assemble mesh velocity*/
         pSgo_rot->Assemble();
-        
+
         /*Assign output pointers:*/
         *pIxz=moi_list[0];
@@ -4412 +4412 @@
         /*Free ressources:*/
         xDelete<IssmDouble>(Sg_old);
-        
+
 }
 /*}}}*/
@@ -4419 +4419 @@
         /*serialized vectors:*/
         IssmDouble* Sg=NULL;
-        
+
         IssmDouble  eartharea=0;
         IssmDouble  eartharea_cpu=0;
 
         int         ns,nsmax;
-        
+
         /*Serialize vectors from previous iteration:*/
         Sg=pSg->ToMPISerial();
@@ -4430 +4430 @@
         /*Go through elements, and add contribution from each element to the deflection vector wg:*/
         ns = elements->Size();
-        
+
         /*First, figure out the area of the ocean, which is needed to compute the eustatic component: */
         for(int i=0;i<ns;i++){
@@ -4455 +4455 @@
                 }
         }
-        
+
         /*One last time: */
         pUp->Assemble();
@@ -4492 +4492 @@
         ISSM_MPI_Reduce (&oceanarea_cpu,&oceanarea,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() );
         ISSM_MPI_Bcast(&oceanarea,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
-        
+
         ISSM_MPI_Reduce (&oceanvalue_cpu,&oceanvalue,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() );
         ISSM_MPI_Bcast(&oceanvalue,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
@@ -4498 +4498 @@
         /*Free ressources:*/
         xDelete<IssmDouble>(Sg_serial);
-        
+
         return oceanvalue/oceanarea;
 }
@@ -4514 +4514 @@
         int*                eplzigzag_counter = NULL;
         int                 eplflip_lock;
-        
+
         HydrologyDCEfficientAnalysis* effanalysis =  new HydrologyDCEfficientAnalysis();
         HydrologyDCInefficientAnalysis* inefanalysis =  new HydrologyDCInefficientAnalysis();
@@ -4521 +4521 @@
         mask=new Vector<IssmDouble>(this->nodes->NumberOfNodes(HydrologyDCEfficientAnalysisEnum));
         recurence=new Vector<IssmDouble>(this->nodes->NumberOfNodes(HydrologyDCEfficientAnalysisEnum));
-        this->parameters->FindParam(&eplzigzag_counter,NULL,EplZigZagCounterEnum); 
-        this->parameters->FindParam(&eplflip_lock,HydrologydcEplflipLockEnum); 
+        this->parameters->FindParam(&eplzigzag_counter,NULL,EplZigZagCounterEnum);
+        this->parameters->FindParam(&eplflip_lock,HydrologydcEplflipLockEnum);
         GetVectorFromInputsx(&old_active,this,HydrologydcMaskEplactiveNodeEnum,NodeSIdEnum);
-        
+
         for (int i=0;i<elements->Size();i++){
                 Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
@@ -4542 +4542 @@
         /*Assemble and serialize*/
         mask->Assemble();
-        serial_mask=mask->ToMPISerial();        
-        
+        serial_mask=mask->ToMPISerial();
+
         xDelete<int>(eplzigzag_counter);
         xDelete<IssmDouble>(serial_rec);
@@ -4585 +4585 @@
         int sum_counter;
         ISSM_MPI_Reduce(&counter,&sum_counter,1,ISSM_MPI_INT,ISSM_MPI_SUM,0,IssmComm::GetComm() );
-        ISSM_MPI_Bcast(&sum_counter,1,ISSM_MPI_INT,0,IssmComm::GetComm());                
+        ISSM_MPI_Bcast(&sum_counter,1,ISSM_MPI_INT,0,IssmComm::GetComm());
         counter=sum_counter;
         *pEplcount = counter;
         if(VerboseSolution()) _printf0_("   Number of active nodes in EPL layer: "<< counter <<"\n");
+
+        /*Update dof indexings*/
+        this->UpdateConstraintsx();
+
+}
+/*}}}*/
+void FemModel::HydrologyIDSupdateDomainx(IssmDouble* pIDScount){ /*{{{*/
+
+        bool                isthermal;
+        Vector<IssmDouble>* mask                                = NULL;
+        IssmDouble*         serial_mask = NULL;
+
+        HydrologyDCInefficientAnalysis* inefanalysis =  new HydrologyDCInefficientAnalysis();
+        parameters->FindParam(&isthermal,TransientIsthermalEnum);
+
+        /*When solving a thermal model we update the thawed nodes*/
+        if(isthermal){
+                /*Step 1: update mask, the mask correspond to thawed nodes (that have a meltingrate)*/
+                mask=new Vector<IssmDouble>(this->nodes->NumberOfNodes(HydrologyDCInefficientAnalysisEnum));
+
+                for (int i=0;i<elements->Size();i++){
+                        Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
+                        inefanalysis->HydrologyIDSGetMask(mask,element);
+                }
+                /*Assemble and serialize*/
+                mask->Assemble();
+                serial_mask=mask->ToMPISerial();
+                delete mask;
+        }
+        /*for other cases we just grab the mask from the initialisation value*/
+        else{
+        GetVectorFromInputsx(&serial_mask,this,HydrologydcMaskThawedNodeEnum,NodeSIdEnum);
+        }
+        /*Update node activation accordingly*/
+        int counter =0;
+        for (int i=0;i<nodes->Size();i++){
+                Node* node=xDynamicCast<Node*>(nodes->GetObjectByOffset(i));
+                if(node->InAnalysis(HydrologyDCInefficientAnalysisEnum)){
+                        if(serial_mask[node->Sid()]==1.){
+                                node->Activate();
+                                if(!node->IsClone()) counter++;
+                        }
+                        else{
+                                node->Deactivate();
+                        }
+                }
+        }
+
+        /*Update Mask*/
+        InputUpdateFromVectorx(this,serial_mask,HydrologydcMaskThawedNodeEnum,NodeSIdEnum);
+        xDelete<IssmDouble>(serial_mask);
+        inefanalysis->ElementizeIdsMask(this);
+
+        int sum_counter;
+        ISSM_MPI_Reduce(&counter,&sum_counter,1,ISSM_MPI_INT,ISSM_MPI_SUM,0,IssmComm::GetComm() );
+        ISSM_MPI_Bcast(&sum_counter,1,ISSM_MPI_INT,0,IssmComm::GetComm());
+        counter=sum_counter;
+        *pIDScount = counter;
+        if(VerboseSolution()) _printf0_("   Number of active nodes in IDS layer: "<< counter <<"\n");
 
         /*Update dof indexings*/
@@ -4631 +4690 @@
         int sum_counter;
         ISSM_MPI_Reduce(&counter,&sum_counter,1,ISSM_MPI_INT,ISSM_MPI_SUM,0,IssmComm::GetComm() );
-        ISSM_MPI_Bcast(&sum_counter,1,ISSM_MPI_INT,0,IssmComm::GetComm());                
+        ISSM_MPI_Bcast(&sum_counter,1,ISSM_MPI_INT,0,IssmComm::GetComm());
         counter=sum_counter;
         *pL2count = counter;
@@ -4716 +4775 @@
 }
 /*}}}*/
-#ifdef _HAVE_JAVASCRIPT_ 
+#ifdef _HAVE_JAVASCRIPT_
 FemModel::FemModel(IssmDouble* buffer, int buffersize, char* toolkits, char* solution, char* modelname,ISSM_MPI_Comm incomm, bool trace){ /*{{{*/
         /*configuration: */
@@ -4731 +4790 @@
         /*From command line arguments, retrieve different filenames needed to create the FemModel: */
         solution_type=StringToEnumx(solution);
-        
+
         /*Create femmodel from input files: */
         profiler->Start(MPROCESSOR);
         this->InitFromBuffers((char*)buffer,buffersize,toolkits, solution_type,trace,NULL);
         profiler->Stop(MPROCESSOR);
-        
+
         /*Save communicator in the parameters dataset: */
         this->parameters->AddObject(new GenericParam<ISSM_MPI_Comm>(incomm,FemModelCommEnum));
@@ -4751 +4810 @@
         size_t* poutputbuffersize;
 
-        
+
         /*Before we delete the profiler, report statistics for this run: */
         profiler->Stop(TOTAL);  //final tagging
@@ -4764 +4823 @@
                                 );
         _printf0_("\n");
-        
+
         /*Before we close the output file, recover the buffer and size:*/
         outputbufferparam = xDynamicCast<GenericParam<char**>*>(this->parameters->FindParamObject(OutputBufferPointerEnum));
@@ -4804 +4863 @@
 
         /*Open output file once for all and add output file descriptor to parameters*/
-        output_fid=open_memstream(&outputbuffer,&outputsize); 
+        output_fid=open_memstream(&outputbuffer,&outputsize);
         if(output_fid==NULL)_error_("could not initialize output stream");
         this->parameters->SetParam(output_fid,OutputFilePointerEnum);
@@ -4845 +4904 @@
                 /*Initialize hmaxvertices with NAN*/
                 hmaxvertices_serial=xNew<IssmDouble>(numberofvertices);
-                for(int i=0;i<numberofvertices;i++) hmaxvertices_serial[i]=NAN; 
+                for(int i=0;i<numberofvertices;i++) hmaxvertices_serial[i]=NAN;
                 /*Fill hmaxvertices*/
                 if(this->amrbamg->groundingline_distance>0)             this->GethmaxVerticesFromZeroLevelSetDistance(hmaxvertices_serial,MaskGroundediceLevelsetEnum);
@@ -4852 +4911 @@
                 if(this->amrbamg->deviatoricerror_threshold>0)  this->GethmaxVerticesFromEstimators(hmaxvertices_serial,DeviatoricStressErrorEstimatorEnum);
         }
-        
+
         if(my_rank==0){
                 this->amrbamg->ExecuteRefinementBamg(vector_serial,hmaxvertices_serial,&newnumberofvertices,&newnumberofelements,&newx,&newy,&newz,&newelementslist);
@@ -4872 +4931 @@
                 newelementslist=xNew<int>(newnumberofelements*this->GetElementsWidth());
         }
-        ISSM_MPI_Bcast(newx,newnumberofvertices,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); 
-        ISSM_MPI_Bcast(newy,newnumberofvertices,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); 
-        ISSM_MPI_Bcast(newz,newnumberofvertices,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); 
-        ISSM_MPI_Bcast(newelementslist,newnumberofelements*this->GetElementsWidth(),ISSM_MPI_INT,0,IssmComm::GetComm());        
+        ISSM_MPI_Bcast(newx,newnumberofvertices,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
+        ISSM_MPI_Bcast(newy,newnumberofvertices,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
+        ISSM_MPI_Bcast(newz,newnumberofvertices,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
+        ISSM_MPI_Bcast(newelementslist,newnumberofelements*this->GetElementsWidth(),ISSM_MPI_INT,0,IssmComm::GetComm());
 
         /*Assign output pointers*/
@@ -4908 +4967 @@
         /*Create bamg data structures for bamg*/
         this->amrbamg = new AmrBamg();
-        
+
         /*Get amr parameters*/
         this->parameters->FindParam(&hmin,AmrHminEnum);
@@ -4932 +4991 @@
 
         /*Re-create original mesh and put it in bamg structure (only cpu 0)*/
-        if(my_rank==0){ 
+        if(my_rank==0){
                 this->amrbamg->Initialize(elements,x,y,numberofvertices,numberofelements);
         }
@@ -4946 +5005 @@
 
         if(!hmaxvertices) _error_("hmaxvertices is NULL!\n");
-        
+
         /*Intermediaries*/
         int numberofvertices                     = this->vertices->NumberOfVertices();
@@ -4953 +5012 @@
 
         switch(levelset_type){
-                case MaskGroundediceLevelsetEnum: 
+                case MaskGroundediceLevelsetEnum:
                         threshold       = this->amrbamg->groundingline_distance;
                         resolution      = this->amrbamg->groundingline_resolution;
@@ -4967 +5026 @@
         this->GetVerticeDistanceToZeroLevelSet(&verticedistance,levelset_type);
         if(!verticedistance) _error_("verticedistance is NULL!\n");
-        
+
         /*Fill hmaxVertices*/
         for(int i=0;i<numberofvertices;i++){
@@ -4981 +5040 @@
 /*}}}*/
 void FemModel::GethmaxVerticesFromEstimators(IssmDouble* hmaxvertices,int errorestimator_type){/*{{{*/
-   
+
         if(!hmaxvertices) _error_("hmaxvertices is NULL!\n");
 
@@ -4989 +5048 @@
         int numberofvertices                                            = this->vertices->NumberOfVertices();
         IssmDouble* maxlength                                   = xNew<IssmDouble>(numberofelements);
-        IssmDouble* error_vertices                              = xNewZeroInit<IssmDouble>(numberofvertices);   
+        IssmDouble* error_vertices                              = xNewZeroInit<IssmDouble>(numberofvertices);
         IssmDouble* error_elements                              = NULL;
         IssmDouble* x                                                           = NULL;
@@ -5002 +5061 @@
         /*Fill variables*/
         switch(errorestimator_type){
-                case ThicknessErrorEstimatorEnum: 
+                case ThicknessErrorEstimatorEnum:
                         threshold               = this->amrbamg->thicknesserror_threshold;
                         groupthreshold  = this->amrbamg->thicknesserror_groupthreshold;
@@ -5027 +5086 @@
         case ThicknessErrorEstimatorEnum:                       this->amrbamg->thicknesserror_maximum   = maxerror;break;
         case DeviatoricStressErrorEstimatorEnum:        this->amrbamg->deviatoricerror_maximum = maxerror;break;
-        }       
+        }
         }
 
         /*Get mesh*/
         this->GetMesh(this->vertices,this->elements,&x,&y,&z,&index);
-        
+
         /*Fill error_vertices (this is the sum of all elements connected to the vertex)*/
         for(int i=0;i<numberofelements;i++){
@@ -5046 +5105 @@
                 error_vertices[v2]+=error_elements[i];
                 error_vertices[v3]+=error_elements[i];
-        }       
+        }
 
         /*Fill hmaxvertices with the criteria*/
@@ -5060 +5119 @@
                         }
                 }
-                /*Now, fill the hmaxvertices if requested*/       
+                /*Now, fill the hmaxvertices if requested*/
                 if(refine){
                         for(int j=0;j<elementswidth;j++){
@@ -5090 +5149 @@
         /*Output*/
         IssmDouble* verticedistance;
-        
+
         /*Intermediaries*/
    int numberofvertices       = this->vertices->NumberOfVertices();
@@ -5102 +5161 @@
         /*Get vertices coordinates*/
         VertexCoordinatesx(&x,&y,&z,this->vertices,false) ;
-        
-        /*Get points which level set is zero (center of elements with zero level set)*/ 
+
+        /*Get points which level set is zero (center of elements with zero level set)*/
         this->GetZeroLevelSetPoints(&levelset_points,numberofpoints,levelset_type);
 
@@ -5112 +5171 @@
                 for(int j=0;j<numberofpoints;j++){
                         distance=sqrt((x[i]-levelset_points[2*j])*(x[i]-levelset_points[2*j])+(y[i]-levelset_points[2*j+1])*(y[i]-levelset_points[2*j+1]));
-                        verticedistance[i]=min(distance,verticedistance[i]);            
-                }
-        }       
+                        verticedistance[i]=min(distance,verticedistance[i]);
+                }
+        }
 
         /*Assign the pointer*/
@@ -5148 +5207 @@
         if(this->amr->groundingline_distance>0)         this->GetElementDistanceToZeroLevelSet(&gl_distance,MaskGroundediceLevelsetEnum);
    if(this->amr->icefront_distance>0)                           this->GetElementDistanceToZeroLevelSet(&if_distance,MaskIceLevelsetEnum);
-   if(this->amr->thicknesserror_threshold>0)            this->ThicknessZZErrorEstimator(&thicknesserror);       
-        if(this->amr->deviatoricerror_threshold>0)      this->ZZErrorEstimator(&deviatoricerror);       
-        
+   if(this->amr->thicknesserror_threshold>0)            this->ThicknessZZErrorEstimator(&thicknesserror);
+        if(this->amr->deviatoricerror_threshold>0)      this->ZZErrorEstimator(&deviatoricerror);
+
         if(my_rank==0){
                 this->amr->ExecuteRefinement(gl_distance,if_distance,deviatoricerror,thicknesserror,
-                                                                                                &newnumberofvertices,&newnumberofelements,&newx,&newy,&newelementslist); 
+                                                                                                &newnumberofvertices,&newnumberofelements,&newx,&newy,&newelementslist);
                 newz=xNewZeroInit<IssmDouble>(newnumberofvertices);
                 if(newnumberofvertices<=0 || newnumberofelements<=0) _error_("Error in the ReMeshNeopz.");
@@ -5167 +5226 @@
                 newelementslist=xNew<int>(newnumberofelements*this->GetElementsWidth());
         }
-        ISSM_MPI_Bcast(newx,newnumberofvertices,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); 
-        ISSM_MPI_Bcast(newy,newnumberofvertices,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); 
-        ISSM_MPI_Bcast(newz,newnumberofvertices,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); 
-        ISSM_MPI_Bcast(newelementslist,newnumberofelements*this->GetElementsWidth(),ISSM_MPI_INT,0,IssmComm::GetComm());        
-
-        /*Assign the pointers*/ 
+        ISSM_MPI_Bcast(newx,newnumberofvertices,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
+        ISSM_MPI_Bcast(newy,newnumberofvertices,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
+        ISSM_MPI_Bcast(newz,newnumberofvertices,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
+        ISSM_MPI_Bcast(newelementslist,newnumberofelements*this->GetElementsWidth(),ISSM_MPI_INT,0,IssmComm::GetComm());
+
+        /*Assign the pointers*/
         (*pnewelementslist)     = newelementslist; //Matlab indexing
         (*pnewx)                                        = newx;
@@ -5188 +5247 @@
 /*}}}*/
 void FemModel::InitializeAdaptiveRefinementNeopz(void){/*{{{*/
-        
+
         /*Define variables*/
         int my_rank                                                                             = IssmComm::GetRank();
@@ -5197 +5256 @@
         IssmDouble* z                                                                   = NULL;
         int* elements                                                                   = NULL;
-        
+
         /*Initialize field as NULL for now*/
         this->amr = NULL;
@@ -5205 +5264 @@
         this->GetMesh(this->vertices,this->elements,&x,&y,&z,&elements);
 
-        /*Create initial mesh (coarse mesh) in neopz data structure*/ 
+        /*Create initial mesh (coarse mesh) in neopz data structure*/
         /*Just CPU #0 should keep AMR object*/
    /*Initialize refinement pattern*/
         this->SetRefPatterns();
         this->amr = new AdaptiveMeshRefinement();
-        this->amr->refinement_type=1;//1 is refpattern; 0 is uniform (faster) 
+        this->amr->refinement_type=1;//1 is refpattern; 0 is uniform (faster)
         /*Get amr parameters*/
         this->parameters->FindParam(&this->amr->level_max,AmrLevelMaxEnum);
@@ -5223 +5282 @@
         this->parameters->FindParam(&this->amr->deviatoricerror_groupthreshold,AmrDeviatoricErrorGroupThresholdEnum);
         this->parameters->FindParam(&this->amr->deviatoricerror_maximum,AmrDeviatoricErrorMaximumEnum);
-        if(my_rank==0){ 
+        if(my_rank==0){
                 this->amr->CreateInitialMesh(numberofvertices,numberofelements,x,y,elements);
         }
@@ -5244 +5303 @@
         /*Output*/
         IssmDouble* elementdistance;
-        
+
         /*Intermediaries*/
    int numberofelements                                                 = this->elements->NumberOfElements();
@@ -5253 +5312 @@
         IssmDouble xc,yc,x1,y1,x2,y2,x3,y3;
         int numberofpoints;
-        
-        /*Get points which level set is zero (center of elements with zero level set, levelset_points is serial)*/      
+
+        /*Get points which level set is zero (center of elements with zero level set, levelset_points is serial)*/
         this->GetZeroLevelSetPoints(&levelset_points,numberofpoints,levelset_type);
 
@@ -5260 +5319 @@
       Element* element=xDynamicCast<Element*>(this->elements->GetObjectByOffset(i));
       int sid = element->Sid();
-                element->GetVerticesCoordinates(&xyz_list); 
+                element->GetVerticesCoordinates(&xyz_list);
                 x1 = xyz_list[3*0+0];y1 = xyz_list[3*0+1];
                 x2 = xyz_list[3*1+0];y2 = xyz_list[3*1+1];
@@ -5270 +5329 @@
                 for(int j=0;j<numberofpoints;j++){
                         distance =sqrt((xc-levelset_points[2*j])*(xc-levelset_points[2*j])+(yc-levelset_points[2*j+1])*(yc-levelset_points[2*j+1]));
-                        mindistance=min(distance,mindistance);          
+                        mindistance=min(distance,mindistance);
                 }
                 velementdistance->SetValue(sid,mindistance,INS_VAL);
                 xDelete<IssmDouble>(xyz_list);
-        }       
+        }
 
    /*Assemble*/

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

@@ -1 +1 @@
 /*
- * FemModel.h: 
+ * FemModel.h:
  */
 
@@ -47 +47 @@
                 Nodes       *nodes;                //one set of nodes
                 Parameters  *parameters;           //one set of parameters, independent of the analysis_type
-                Results     *results;              //results that cannot be fit into the elements 
+                Results     *results;              //results that cannot be fit into the elements
                 Vertices    *vertices;             //one set of vertices
 
@@ -80 +80 @@
                 void Solve(void);
 
-                /*Modules*/ 
+                /*Modules*/
                 void BalancethicknessMisfitx(IssmDouble* pV);
                 void CalvingRateVonmisesx();
@@ -135 +135 @@
                 #endif
                 #ifdef _HAVE_ESA_
-                void EsaGeodetic2D(Vector<IssmDouble>* pUp, Vector<IssmDouble>* pNorth, Vector<IssmDouble>* pEast, Vector<IssmDouble>* pX, Vector<IssmDouble>* pY, IssmDouble* xx, IssmDouble* yy); 
-                void EsaGeodetic3D(Vector<IssmDouble>* pUp, Vector<IssmDouble>* pNorth, Vector<IssmDouble>* pEast, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz); 
+                void EsaGeodetic2D(Vector<IssmDouble>* pUp, Vector<IssmDouble>* pNorth, Vector<IssmDouble>* pEast, Vector<IssmDouble>* pX, Vector<IssmDouble>* pY, IssmDouble* xx, IssmDouble* yy);
+                void EsaGeodetic3D(Vector<IssmDouble>* pUp, Vector<IssmDouble>* pNorth, Vector<IssmDouble>* pEast, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz);
                 #endif
                 #ifdef _HAVE_SEALEVELRISE_
@@ -142 +142 @@
                 void SealevelriseNonEustatic(Vector<IssmDouble>* pSgo, Vector<IssmDouble>* pSg_old, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius,bool verboseconvolution);
                 void SealevelriseRotationalFeedback(Vector<IssmDouble>* pSgo_rot, Vector<IssmDouble>* pSg_old, IssmDouble* pIxz, IssmDouble* pIyz, IssmDouble* pIzz, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius);
-                void SealevelriseGeodetic(Vector<IssmDouble>* pUp, Vector<IssmDouble>* pNorth, Vector<IssmDouble>* pEast, Vector<IssmDouble>* pSg_old, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz); 
+                void SealevelriseGeodetic(Vector<IssmDouble>* pUp, Vector<IssmDouble>* pNorth, Vector<IssmDouble>* pEast, Vector<IssmDouble>* pSg_old, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz);
                 IssmDouble SealevelriseOceanAverage(Vector<IssmDouble>* Sg);
                 #endif
                 void HydrologyEPLupdateDomainx(IssmDouble* pEplcount);
+                void HydrologyIDSupdateDomainx(IssmDouble* pIDScount);
                 void TimeAdaptx(IssmDouble* pdt);
                 void UpdateConstraintsExtrudeFromBasex();

issm/trunk-jpl/src/c/cores/hydrology_core.cpp

@@ -13 +13 @@
 
         /*intermediary*/
-        int  hydrology_model;
-        int  solution_type;
-        int  numoutputs=0;
-        bool save_results;
-        bool modify_loads=true;
-        char **requested_outputs = NULL;
+        int                                                     hydrology_model;
+        int                                                     solution_type;
+        int                                                     numoutputs                              =       0;
+        bool                                            save_results;
+        bool                                            modify_loads                    =       true;
+        char                                    **requested_outputs = NULL;
+        IssmDouble                      ThawedNodes;
 
         /*first recover parameters common to all solutions*/
@@ -45 +46 @@
         else if (hydrology_model==HydrologydcEnum){
                 /*intermediary: */
-                bool       isefficientlayer;
-                int        step,hydroslices;
-                IssmDouble time,init_time,hydrotime,yts;
-                IssmDouble dt,hydrodt;
+                bool                            isefficientlayer;
+                bool                            isthermal;
+                int                                     step,hydroslices;
+                IssmDouble      time,init_time,hydrotime,yts;
+                IssmDouble      dt,hydrodt;
 
                 femmodel->parameters->FindParam(&isefficientlayer,HydrologydcIsefficientlayerEnum);
@@ -57 +59 @@
                 femmodel->parameters->FindParam(&yts,ConstantsYtsEnum);
 
-                init_time = time-dt; //getting the time back to the start of the timestep
-                hydrotime=init_time;
-                hydrodt=dt/hydroslices; //computing hydro dt from dt and a divider
-                femmodel->parameters->AddObject(new DoubleParam(HydrologydtEnum,hydrodt));
-                if(hydroslices>1){
-                        /*define which variable needs to be averaged on the sub-timestep and initialize as needed*/
-                        if (isefficientlayer){
-                                int inputtostack[4]={EffectivePressureHydrostepEnum,SedimentHeadHydrostepEnum,EplHeadHydrostepEnum,HydrologydcEplThicknessHydrostepEnum};
-                                int stackedinput[4]={EffectivePressureStackedEnum,SedimentHeadStackedEnum,EplHeadStackedEnum,HydrologydcEplThicknessStackedEnum};
-                                int averagedinput[4]={EffectivePressureEnum,SedimentHeadEnum,EplHeadEnum,HydrologydcEplThicknessEnum};
-                                femmodel->InitTransientOutputx(&stackedinput[0],4);
-                                while(hydrotime<time-(yts*DBL_EPSILON)){ //loop on hydro dts
-                                        hydrotime+=hydrodt;
-                                        /*save preceding timestep*/
-                                        InputDuplicatex(femmodel,SedimentHeadHydrostepEnum,SedimentHeadOldEnum);
+                femmodel->parameters->FindParam(&isthermal,TransientIsthermalEnum);
+                /*first we exclude frozen nodes of the solved nodes*/
+                femmodel->SetCurrentConfiguration(HydrologyDCInefficientAnalysisEnum);
+                femmodel->HydrologyIDSupdateDomainx(&ThawedNodes);
+                if(ThawedNodes>0){
+                        init_time = time-dt; //getting the time back to the start of the timestep
+                        hydrotime=init_time;
+                        hydrodt=dt/hydroslices; //computing hydro dt from dt and a divider
+                        femmodel->parameters->AddObject(new DoubleParam(HydrologydtEnum,hydrodt));
+                        if(hydroslices>1){
+                                /*define which variable needs to be averaged on the sub-timestep and initialize as needed*/
+                                if (isefficientlayer){
+                                        int inputtostack[4]={EffectivePressureHydrostepEnum,SedimentHeadHydrostepEnum,EplHeadHydrostepEnum,HydrologydcEplThicknessHydrostepEnum};
+                                        int stackedinput[4]={EffectivePressureStackedEnum,SedimentHeadStackedEnum,EplHeadStackedEnum,HydrologydcEplThicknessStackedEnum};
+                                        int averagedinput[4]={EffectivePressureEnum,SedimentHeadEnum,EplHeadEnum,HydrologydcEplThicknessEnum};
+                                        femmodel->InitTransientOutputx(&stackedinput[0],4);
+                                        while(hydrotime<time-(yts*DBL_EPSILON)){ //loop on hydro dts
+                                                hydrotime+=hydrodt;
+                                                /*save preceding timestep*/
+                                                InputDuplicatex(femmodel,SedimentHeadHydrostepEnum,SedimentHeadOldEnum);
+                                                InputDuplicatex(femmodel,EplHeadHydrostepEnum,EplHeadOldEnum);
+                                                InputDuplicatex(femmodel,HydrologydcEplThicknessHydrostepEnum,HydrologydcEplThicknessOldEnum);
+                                                /*Proceed now to heads computations*/
+                                                solutionsequence_hydro_nonlinear(femmodel);
+                                                /*If we have a sub-timestep we stack the variables here*/
+                                                femmodel->StackTransientOutputx(&inputtostack[0],&stackedinput[0],hydrotime,4);
+                                        }
+                                        femmodel->AverageTransientOutputx(&stackedinput[0],&averagedinput[0],init_time,4);
+                                }
+                                else{
+                                        int inputtostack[2]={EffectivePressureHydrostepEnum,SedimentHeadHydrostepEnum};
+                                        int stackedinput[2]={EffectivePressureStackedEnum,SedimentHeadStackedEnum};
+                                        int averagedinput[2]={EffectivePressureEnum,SedimentHeadEnum};
+                                        femmodel->InitTransientOutputx(&stackedinput[0],2);
+                                        while(hydrotime<time-(yts*DBL_EPSILON)){ //loop on hydro dts
+                                                hydrotime+=hydrodt;
+                                                /*save preceding timestep*/
+                                                InputDuplicatex(femmodel,SedimentHeadHydrostepEnum,SedimentHeadOldEnum);
+                                                /*Proceed now to heads computations*/
+                                                solutionsequence_hydro_nonlinear(femmodel);
+                                                /*If we have a sub-timestep we stack the variables here*/
+                                                femmodel->StackTransientOutputx(&inputtostack[0],&stackedinput[0],hydrotime,2);
+                                        }
+                                        femmodel->AverageTransientOutputx(&stackedinput[0],&averagedinput[0],init_time,2);
+                                }
+                        }
+                        else{
+                                InputDuplicatex(femmodel,SedimentHeadHydrostepEnum,SedimentHeadOldEnum);
+                                if (isefficientlayer){
                                         InputDuplicatex(femmodel,EplHeadHydrostepEnum,EplHeadOldEnum);
                                         InputDuplicatex(femmodel,HydrologydcEplThicknessHydrostepEnum,HydrologydcEplThicknessOldEnum);
-                                        /*Proceed now to heads computations*/
-                                        solutionsequence_hydro_nonlinear(femmodel);
-                                        /*If we have a sub-timestep we stack the variables here*/
-                                        femmodel->StackTransientOutputx(&inputtostack[0],&stackedinput[0],hydrotime,4);
                                 }
-                                femmodel->AverageTransientOutputx(&stackedinput[0],&averagedinput[0],init_time,4);
+                                /*Proceed now to heads computations*/
+                                solutionsequence_hydro_nonlinear(femmodel);
                         }
-                        else{
-                                int inputtostack[2]={EffectivePressureHydrostepEnum,SedimentHeadHydrostepEnum};
-                                int stackedinput[2]={EffectivePressureStackedEnum,SedimentHeadStackedEnum};
-                                int averagedinput[2]={EffectivePressureEnum,SedimentHeadEnum};
-                                femmodel->InitTransientOutputx(&stackedinput[0],2);
-                                while(hydrotime<time-(yts*DBL_EPSILON)){ //loop on hydro dts
-                                        hydrotime+=hydrodt;
-                                        /*save preceding timestep*/
-                                        InputDuplicatex(femmodel,SedimentHeadHydrostepEnum,SedimentHeadOldEnum);
-                                        /*Proceed now to heads computations*/
-                                        solutionsequence_hydro_nonlinear(femmodel);
-                                        /*If we have a sub-timestep we stack the variables here*/
-                                        femmodel->StackTransientOutputx(&inputtostack[0],&stackedinput[0],hydrotime,2);
-                                }
-                                femmodel->AverageTransientOutputx(&stackedinput[0],&averagedinput[0],init_time,2);
-                        }
-                }
-                else{
-                        InputDuplicatex(femmodel,SedimentHeadHydrostepEnum,SedimentHeadOldEnum);
-                        if (isefficientlayer){
-                                InputDuplicatex(femmodel,EplHeadHydrostepEnum,EplHeadOldEnum);
-                                InputDuplicatex(femmodel,HydrologydcEplThicknessHydrostepEnum,HydrologydcEplThicknessOldEnum);
-                        }
-                        /*Proceed now to heads computations*/
-                        solutionsequence_hydro_nonlinear(femmodel);
                 }
         }
         else if (hydrology_model==HydrologysommersEnum){
                 femmodel->SetCurrentConfiguration(HydrologySommersAnalysisEnum);
-      InputDuplicatex(femmodel,HydrologyHeadEnum,HydrologyHeadOldEnum);
+                InputDuplicatex(femmodel,HydrologyHeadEnum,HydrologyHeadOldEnum);
                 solutionsequence_shakti_nonlinear(femmodel);
                 if(VerboseSolution()) _printf0_("   updating gap height\n");
@@ -122 +130 @@
         if(save_results){
                 if(VerboseSolution()) _printf0_("   saving results \n");
-                femmodel->RequestedOutputsx(&femmodel->results,requested_outputs,numoutputs);
+                if(hydrology_model==HydrologydcEnum && ThawedNodes==0){
+                        if(VerboseSolution()) _printf0_("   No thawed node hydro is skiped \n");}
+                else{
+                        femmodel->RequestedOutputsx(&femmodel->results,requested_outputs,numoutputs);
+                }
         }
         /*Free ressources:*/

issm/trunk-jpl/src/c/cores/transient_core.cpp

@@ -1 +1 @@
 /*!\file: transient_3d_core.cpp
- * \brief: core of the transient_3d solution 
- */ 
+ * \brief: core of the transient_3d solution
+ */
 
 #ifdef HAVE_CONFIG_H
@@ -142 +142 @@
                         InterpFromMeshToMesh2dx(&base_oceangrid,index_ice,lon_ice,lat_ice,ngrids_ice,nels_ice,
                                                         icebase,ngrids_ice,1,
-                                                        oceangridx,oceangridy,ngrids_ocean,options);   
+                                                        oceangridx,oceangridy,ngrids_ocean,options);
                         delete options;
 
@@ -239 +239 @@
                                 int         ngrids_ice=femmodel->vertices->NumberOfVertices();
                                 int         nels_ice=femmodel->elements->NumberOfElements();
-                                
+
                                 /*Recover mesh and inputs needed*/
                                 femmodel->parameters->FindParam(&rho_ice,MaterialsRhoIceEnum);
@@ -253 +253 @@
                                 InterpFromMeshToMesh2dx(&base_oceangrid,index_ice,lon_ice,lat_ice,ngrids_ice,nels_ice,
                                                         icebase,ngrids_ice,1,
-                                                        oceangridx,oceangridy,ngrids_ocean,NULL);   
-                                
+                                                        oceangridx,oceangridy,ngrids_ocean,NULL);
+
                                 /*Send and receive data*/
                                 ISSM_MPI_Send(&time,1,ISSM_MPI_DOUBLE,0,10001001,tomitgcmcomm);
@@ -267 +267 @@
                                 InterpFromMeshToMesh2dx(&melt_mesh,index_ocean,oceangridx,oceangridy,ngrids_ocean,nels_ocean,
                                                         oceanmelt,ngrids_ocean,1,
-                                                        lon_ice,lat_ice,ngrids_ice,NULL);   
+                                                        lon_ice,lat_ice,ngrids_ice,NULL);
 
                                 for(int i=0;i<ngrids_ice;i++) melt_mesh[i]=-melt_mesh[i]/rho_ice; //heat flux provided by ocean is in kg/m^2/s
@@ -293 +293 @@
                 /*}}}*/
 
-                if(isthermal && domaintype==Domain3DEnum){ 
+                if(isthermal && domaintype==Domain3DEnum){
                         if(issmb){
                                 bool isenthalpy;
@@ -327 +327 @@
                         femmodel->UpdateVertexPositionsx();
                 }
-                
+
                 if(isgroundingline){
                         if(VerboseSolution()) _printf0_("   computing new grounding line position\n");
@@ -338 +338 @@
                         femmodel->parameters->SetParam(SurfaceEnum,InputToExtrudeEnum);
                         extrudefrombase_core(femmodel);
-                                
+
                         if(save_results){
                                 int outputs[3] = {SurfaceEnum,BaseEnum,MaskGroundediceLevelsetEnum};
@@ -356 +356 @@
                 /*esa: */
                 if(isesa) esa_core(femmodel);
-                
+
                 /*Sea level rise: */
                 if(isslr || iscoupler) sealevelrise_core(femmodel);
@@ -367 +367 @@
                         femmodel->RequestedOutputsx(&femmodel->results,&outputs[0],1,save_results);
                 }
-                
+
                 if(save_results){
                         if(VerboseSolution()) _printf0_("   saving temporary results\n");
@@ -402 +402 @@
                         }
         }
-        
+
         if(!iscontrol || !isautodiff) femmodel->RequestedDependentsx();
         if(iscontrol && isautodiff) femmodel->parameters->SetParam(dependent_objects,AutodiffDependentObjectsEnum);
 
-        /*Free ressources:*/    
+        /*Free ressources:*/
         if(numoutputs){for(int i=0;i<numoutputs;i++){xDelete<char>(requested_outputs[i]);} xDelete<char*>(requested_outputs);}
 }

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

@@ -441 +441 @@
         HydrologydcMaskEplactiveEltEnum,
         HydrologydcMaskEplactiveNodeEnum,
+        HydrologydcMaskThawedEltEnum,
+        HydrologydcMaskThawedNodeEnum,
         HydrologydcSedimentTransmitivityEnum,
         HydrologyEnglacialInputEnum,

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

@@ -447 +447 @@
                 case HydrologydcMaskEplactiveEltEnum : return "HydrologydcMaskEplactiveElt";
                 case HydrologydcMaskEplactiveNodeEnum : return "HydrologydcMaskEplactiveNode";
+                case HydrologydcMaskThawedEltEnum : return "HydrologydcMaskThawedElt";
+                case HydrologydcMaskThawedNodeEnum : return "HydrologydcMaskThawedNode";
                 case HydrologydcSedimentTransmitivityEnum : return "HydrologydcSedimentTransmitivity";
                 case HydrologyEnglacialInputEnum : return "HydrologyEnglacialInput";

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

@@ -456 +456 @@
               else if (strcmp(name,"HydrologydcMaskEplactiveElt")==0) return HydrologydcMaskEplactiveEltEnum;
               else if (strcmp(name,"HydrologydcMaskEplactiveNode")==0) return HydrologydcMaskEplactiveNodeEnum;
+              else if (strcmp(name,"HydrologydcMaskThawedElt")==0) return HydrologydcMaskThawedEltEnum;
+              else if (strcmp(name,"HydrologydcMaskThawedNode")==0) return HydrologydcMaskThawedNodeEnum;
               else if (strcmp(name,"HydrologydcSedimentTransmitivity")==0) return HydrologydcSedimentTransmitivityEnum;
               else if (strcmp(name,"HydrologyEnglacialInput")==0) return HydrologyEnglacialInputEnum;
@@ -504 +506 @@
               else if (strcmp(name,"Pressure")==0) return PressureEnum;
               else if (strcmp(name,"RheologyBAbsGradient")==0) return RheologyBAbsGradientEnum;
-              else if (strcmp(name,"RheologyBbarAbsGradient")==0) return RheologyBbarAbsGradientEnum;
-              else if (strcmp(name,"Sealevel")==0) return SealevelEnum;
          else stage=5;
    }
    if(stage==5){
-              if (strcmp(name,"SealevelriseDeltathickness")==0) return SealevelriseDeltathicknessEnum;
+              if (strcmp(name,"RheologyBbarAbsGradient")==0) return RheologyBbarAbsGradientEnum;
+              else if (strcmp(name,"Sealevel")==0) return SealevelEnum;
+              else if (strcmp(name,"SealevelriseDeltathickness")==0) return SealevelriseDeltathicknessEnum;
               else if (strcmp(name,"SedimentHeadHydrostep")==0) return SedimentHeadHydrostepEnum;
               else if (strcmp(name,"SedimentHeadOld")==0) return SedimentHeadOldEnum;
@@ -627 +629 @@
               else if (strcmp(name,"VzMesh")==0) return VzMeshEnum;
               else if (strcmp(name,"VzSSA")==0) return VzSSAEnum;
-              else if (strcmp(name,"Watercolumn")==0) return WatercolumnEnum;
-              else if (strcmp(name,"WaterColumnOld")==0) return WaterColumnOldEnum;
          else stage=6;
    }
    if(stage==6){
-              if (strcmp(name,"WaterfractionDrainage")==0) return WaterfractionDrainageEnum;
+              if (strcmp(name,"Watercolumn")==0) return WatercolumnEnum;
+              else if (strcmp(name,"WaterColumnOld")==0) return WaterColumnOldEnum;
+              else if (strcmp(name,"WaterfractionDrainage")==0) return WaterfractionDrainageEnum;
               else if (strcmp(name,"WaterfractionDrainageIntegrated")==0) return WaterfractionDrainageIntegratedEnum;
               else if (strcmp(name,"Waterfraction")==0) return WaterfractionEnum;
@@ -750 +752 @@
               else if (strcmp(name,"FSpressure")==0) return FSpressureEnum;
               else if (strcmp(name,"FSSolver")==0) return FSSolverEnum;
-              else if (strcmp(name,"FSvelocity")==0) return FSvelocityEnum;
-              else if (strcmp(name,"GaussPenta")==0) return GaussPentaEnum;
          else stage=7;
    }
    if(stage==7){
-              if (strcmp(name,"GaussSeg")==0) return GaussSegEnum;
+              if (strcmp(name,"FSvelocity")==0) return FSvelocityEnum;
+              else if (strcmp(name,"GaussPenta")==0) return GaussPentaEnum;
+              else if (strcmp(name,"GaussSeg")==0) return GaussSegEnum;
               else if (strcmp(name,"GaussTetra")==0) return GaussTetraEnum;
               else if (strcmp(name,"GaussTria")==0) return GaussTriaEnum;
@@ -873 +875 @@
               else if (strcmp(name,"MaxAbsVz")==0) return MaxAbsVzEnum;
               else if (strcmp(name,"MaxDivergence")==0) return MaxDivergenceEnum;
-              else if (strcmp(name,"MaxVel")==0) return MaxVelEnum;
-              else if (strcmp(name,"MaxVx")==0) return MaxVxEnum;
          else stage=8;
    }
    if(stage==8){
-              if (strcmp(name,"MaxVy")==0) return MaxVyEnum;
+              if (strcmp(name,"MaxVel")==0) return MaxVelEnum;
+              else if (strcmp(name,"MaxVx")==0) return MaxVxEnum;
+              else if (strcmp(name,"MaxVy")==0) return MaxVyEnum;
               else if (strcmp(name,"MaxVz")==0) return MaxVzEnum;
               else if (strcmp(name,"Melange")==0) return MelangeEnum;
@@ -996 +998 @@
               else if (strcmp(name,"Outputdefinition88")==0) return Outputdefinition88Enum;
               else if (strcmp(name,"Outputdefinition89")==0) return Outputdefinition89Enum;
-              else if (strcmp(name,"Outputdefinition8")==0) return Outputdefinition8Enum;
-              else if (strcmp(name,"Outputdefinition90")==0) return Outputdefinition90Enum;
          else stage=9;
    }
    if(stage==9){
-              if (strcmp(name,"Outputdefinition91")==0) return Outputdefinition91Enum;
+              if (strcmp(name,"Outputdefinition8")==0) return Outputdefinition8Enum;
+              else if (strcmp(name,"Outputdefinition90")==0) return Outputdefinition90Enum;
+              else if (strcmp(name,"Outputdefinition91")==0) return Outputdefinition91Enum;
               else if (strcmp(name,"Outputdefinition92")==0) return Outputdefinition92Enum;
               else if (strcmp(name,"Outputdefinition93")==0) return Outputdefinition93Enum;
@@ -1119 +1121 @@
               else if (strcmp(name,"TotalSmbScaled")==0) return TotalSmbScaledEnum;
               else if (strcmp(name,"TransientArrayParam")==0) return TransientArrayParamEnum;
-              else if (strcmp(name,"TransientInput")==0) return TransientInputEnum;
-              else if (strcmp(name,"TransientParam")==0) return TransientParamEnum;
          else stage=10;
    }
    if(stage==10){
-              if (strcmp(name,"TransientSolution")==0) return TransientSolutionEnum;
+              if (strcmp(name,"TransientInput")==0) return TransientInputEnum;
+              else if (strcmp(name,"TransientParam")==0) return TransientParamEnum;
+              else if (strcmp(name,"TransientSolution")==0) return TransientSolutionEnum;
               else if (strcmp(name,"Tria")==0) return TriaEnum;
               else if (strcmp(name,"TriaInput")==0) return TriaInputEnum;

issm/trunk-jpl/src/c/solutionsequences/solutionsequence_hydro_nonlinear.cpp

@@ -1 +1 @@
 /*!\file: solutionsequence_hydro_nonlinear.cpp
- * \brief: core of the hydro solution 
- */ 
+ * \brief: core of the hydro solution
+ */
 
 #include "../toolkits/toolkits.h"
@@ -12 +12 @@
 void solutionsequence_hydro_nonlinear(FemModel* femmodel){
         /*solution : */
-        Vector<IssmDouble>* ug_sed=NULL; 
-        Vector<IssmDouble>* uf_sed=NULL; 
-        Vector<IssmDouble>* uf_sed_sub_iter=NULL; 
+        Vector<IssmDouble>* ug_sed=NULL;
+        Vector<IssmDouble>* uf_sed=NULL;
+        Vector<IssmDouble>* uf_sed_sub_iter=NULL;
         Vector<IssmDouble>* ug_sed_main_iter=NULL;
 
-        Vector<IssmDouble>* ug_epl=NULL; 
+        Vector<IssmDouble>* ug_epl=NULL;
         Vector<IssmDouble>* uf_epl=NULL;
-        Vector<IssmDouble>* uf_epl_sub_iter=NULL; 
+        Vector<IssmDouble>* uf_epl_sub_iter=NULL;
         Vector<IssmDouble>* ug_epl_sub_iter=NULL;
         Vector<IssmDouble>* ug_epl_main_iter=NULL;
 
-        Vector<IssmDouble>* ys=NULL; 
+        Vector<IssmDouble>* ys=NULL;
         Vector<IssmDouble>* dug=NULL;
         Vector<IssmDouble>* duf=NULL;
@@ -35 +35 @@
         HydrologyDCInefficientAnalysis* inefanalysis = NULL;
         HydrologyDCEfficientAnalysis* effanalysis = NULL;
-        
+
         bool       sedconverged,eplconverged,hydroconverged;
         bool       isefficientlayer;
@@ -58 +58 @@
 
         /*Retrieve inputs as the initial state for the non linear iteration*/
-        GetSolutionFromInputsx(&ug_sed,femmodel);       
+        GetSolutionFromInputsx(&ug_sed,femmodel);
         Reducevectorgtofx(&uf_sed, ug_sed, femmodel->nodes,femmodel->parameters);
+        /*Initialize the thawed element mask*/
         if(isefficientlayer) {
                 inefanalysis = new HydrologyDCInefficientAnalysis();
@@ -87 +88 @@
                 InputUpdateFromConstantx(femmodel,false,ConvergedEnum);
                 femmodel->UpdateConstraintsx();
-                
+
                 /*Reset constraint on the ZigZag Lock*/
                 ResetConstraintsx(femmodel);
@@ -126 +127 @@
                                 if(sedconverged)break;
                         }
-                
+
                         /*}}}*//*End of the sediment penalization loop*/
                         sedconverged=false;
-                        
+
                         /*Checking convegence on the value of the sediment head*/
                         duf=uf_sed_sub_iter->Duplicate();_assert_(duf);
@@ -184 +185 @@
                                 inefanalysis->ElementizeEplMask(femmodel);
                                 /*}}}*/
-                                        
+
                                 if(VerboseSolution()) _printf0_("Building EPL Matrix...\n");
                                 SystemMatricesx(&Kff,&Kfs,&pf,&df,NULL,femmodel);
@@ -197 +198 @@
                                 uf_epl_sub_iter=uf_epl->Duplicate();
                                 uf_epl->Copy(uf_epl_sub_iter);
-                                delete ug_epl; 
+                                delete ug_epl;
                                 Mergesolutionfromftogx(&ug_epl,uf_epl,ys,femmodel->nodes,femmodel->parameters); delete ys;
                                 InputUpdateFromSolutionx(femmodel,ug_epl);
                                 ConstraintsStatex(&constraints_converged,&num_unstable_constraints,femmodel);
-                                                
+
                                 dug=ug_epl_sub_iter->Duplicate();_assert_(dug);
                                 ug_epl_sub_iter->Copy(dug);
@@ -218 +219 @@
                                 /* if(ThickCount<L2Count)eplconverged=true; */
                                 eplcount++;
-                                
+
                                 delete ug_epl_sub_iter;
                                 if(eplconverged){
@@ -235 +236 @@
                         //compute norm(du)/norm(u)
                         dug=ug_sed_main_iter->Duplicate(); _assert_(dug);
-                        ug_sed_main_iter->Copy(dug);    
+                        ug_sed_main_iter->Copy(dug);
                         dug->AYPX(ug_sed,-1.0);
-                        ndu_sed=dug->Norm(NORM_TWO); 
+                        ndu_sed=dug->Norm(NORM_TWO);
                         delete dug;
                         nu_sed=ug_sed_main_iter->Norm(NORM_TWO);
@@ -243 +244 @@
                         if (xIsNan<IssmDouble>(ndu_sed) || xIsNan<IssmDouble>(nu_sed)) _error_("Sed convergence criterion is NaN!");
                         if (ndu_sed==0.0 && nu_sed==0.0) nu_sed=1.0e-6; /*Hacking the case where the Sediment is used but empty*/
-                        dug=ug_epl_main_iter->Duplicate();_assert_(dug); 
-                        ug_epl_main_iter->Copy(dug); 
+                        dug=ug_epl_main_iter->Duplicate();_assert_(dug);
+                        ug_epl_main_iter->Copy(dug);
                         dug->AYPX(ug_epl,-1.0);
-                        ndu_epl=dug->Norm(NORM_TWO); 
+                        ndu_epl=dug->Norm(NORM_TWO);
                         delete dug;
                         nu_epl=ug_epl_main_iter->Norm(NORM_TWO);
@@ -257 +258 @@
                                         hydroconverged=true;
                                 }
-                                else{ 
+                                else{
                                         if(VerboseConvergence()) _printf0_(setw(50) << left << "   Sediment Convergence criterion:" << ndu_sed/nu_sed*100 << "%, aiming lower than " << eps_hyd*100 << " %\n");
                                         if(VerboseConvergence()) _printf0_(setw(50) << left << "   EPL Convergence criterion:" << ndu_epl/nu_epl*100 << "%, aiming lower than " << eps_hyd*100 << " %\n");