Changeset 22955 for issm/trunk-jpl/src/c/classes/FemModel.cpp

Timestamp:

07/16/18 15:39:26 (7 years ago)

Author:

Eric.Larour

Message:

CHG: new reorganization of the slr capabilities.

File:

• issm/trunk-jpl/src/c/classes/FemModel.cpp (modified) (141 diffs)

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

@@ -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;
@@ -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();
@@ -4226 +4226 @@
 #endif
 #ifdef _HAVE_SEALEVELRISE_
-void FemModel::SealevelriseEustatic(Vector<IssmDouble>* pSgi, IssmDouble* peustatic, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius) { /*{{{*/
+void FemModel::SealevelriseEustatic(Vector<IssmDouble>* pRSLgi, IssmDouble* peustatic, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius,int loop) { /*{{{*/
 
         /*serialized vectors:*/
@@ -4262 +4262 @@
                 if(i<ns){
 
-                        if(VerboseConvergence())if(i%100==0)_printf0_("\r" << "      convolution progress: " << (double)i/(double)ns*100 << "%  ");
+                        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);
+                        element->SealevelriseEustatic(pRSLgi,&eustatic_cpu_e,latitude,longitude,radius,oceanarea,eartharea);
                         eustatic_cpu+=eustatic_cpu_e;
                 }
-                if(i%100==0)pSgi->Assemble();
+                if(i%loop==0)pRSLgi->Assemble();
         }
         if(VerboseConvergence())_printf0_("\n");
-
+                
         /*One last time: */
-        pSgi->Assemble();
+        pRSLgi->Assemble();
 
         /*Sum all eustatic components from all cpus:*/
@@ -4285 +4285 @@
 }
 /*}}}*/
-void FemModel::SealevelriseNonEustatic(Vector<IssmDouble>* pSgo, Vector<IssmDouble>* pSg_old, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius, bool verboseconvolution){/*{{{*/
+void FemModel::SealevelriseNonEustatic(Vector<IssmDouble>* pRSLgo, Vector<IssmDouble>* pRSLg_old, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius, bool verboseconvolution,int loop){/*{{{*/
 
         /*serialized vectors:*/
-        IssmDouble* Sg_old=NULL;
-
+        IssmDouble* RSLg_old=NULL;
+        
         IssmDouble  eartharea=0;
         IssmDouble  eartharea_cpu=0;
 
         int         ns,nsmax;
-
+        
         /*Serialize vectors from previous iteration:*/
-        Sg_old=pSg_old->ToMPISerial();
+        RSLg_old=pRSLg_old->ToMPISerial();
 
         /*Go through elements, and add contribution from each element to the deflection vector wg:*/
@@ -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());
@@ -4317 +4317 @@
         for(int i=0;i<nsmax;i++){
                 if(i<ns){
-                        if(verboseconvolution)if(VerboseConvergence())if(i%100==0)_printf_("\r" << "      convolution progress: " << (double)i/(double)ns*100 << "%   ");
+                        if(verboseconvolution)if(VerboseConvergence())if(i%100==0)_printf0_("\r" << "              convolution progress: " << (double)i/(double)ns*100 << "%   ");
                         Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
-                        element->SealevelriseNonEustatic(pSgo,Sg_old,latitude,longitude,radius,eartharea);
-                }
-                if(i%100==0)pSgo->Assemble();
-        }
-        if(verboseconvolution)if(VerboseConvergence())_printf_("\n");
-
+                        element->SealevelriseNonEustatic(pRSLgo,RSLg_old,latitude,longitude,radius,eartharea);
+                }
+                if(i%loop==0)pRSLgo->Assemble();
+        }
+        if(verboseconvolution)if(VerboseConvergence())_printf0_("\n");
+        
         /*Free ressources:*/
-        xDelete<IssmDouble>(Sg_old);
-}
-/*}}}*/
-void FemModel::SealevelriseRotationalFeedback(Vector<IssmDouble>* pSgo_rot, Vector<IssmDouble>* pSg_old, IssmDouble* pIxz, IssmDouble* pIyz, IssmDouble* pIzz, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius){/*{{{*/
+        xDelete<IssmDouble>(RSLg_old);
+}
+/*}}}*/
+void FemModel::SealevelriseRotationalFeedback(Vector<IssmDouble>* pRSLgo_rot, Vector<IssmDouble>* pRSLg_old, IssmDouble* pIxz, IssmDouble* pIyz, IssmDouble* pIzz, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius){/*{{{*/
 
         /*serialized vectors:*/
-        IssmDouble* Sg_old=NULL;
+        IssmDouble* RSLg_old=NULL;
         IssmDouble  eartharea=0;
         IssmDouble  eartharea_cpu=0;
@@ -4341 +4341 @@
 
         /*Serialize vectors from previous iteration:*/
-        Sg_old=pSg_old->ToMPISerial();
+        RSLg_old=pRSLg_old->ToMPISerial();
 
         /*First, figure out the area of the ocean, which is needed to compute the eustatic component: */
@@ -4355 +4355 @@
         for(int i=0;i<elements->Size();i++){
                 Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
-                element->SealevelriseMomentOfInertia(&moi_list[0],Sg_old,eartharea);
+                element->SealevelriseMomentOfInertia(&moi_list[0],RSLg_old,eartharea);
                 moi_list_cpu[0] += moi_list[0];
                 moi_list_cpu[1] += moi_list[1];
@@ -4399 +4399 @@
                                                 (-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
+                pRSLgo_rot->SetValue(sid,value,INS_VAL); //INS_VAL ensures that you don't add several times
         }
 
         /*Assemble mesh velocity*/
-        pSgo_rot->Assemble();
+        pRSLgo_rot->Assemble();
 
         /*Assign output pointers:*/
-        *pIxz=moi_list[0];
-        *pIyz=moi_list[1];
-        *pIzz=moi_list[2];
+        if(pIxz)*pIxz=moi_list[0];
+        if(pIyz)*pIyz=moi_list[1];
+        if(pIzz)*pIzz=moi_list[2];
 
         /*Free ressources:*/
-        xDelete<IssmDouble>(Sg_old);
-
-}
-/*}}}*/
-void FemModel::SealevelriseGeodetic(Vector<IssmDouble>* pUp, Vector<IssmDouble>* pNorth, Vector<IssmDouble>* pEast, Vector<IssmDouble>* pSg, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz){/*{{{*/
+        xDelete<IssmDouble>(RSLg_old);
+        
+        
+}
+/*}}}*/
+void FemModel::SealevelriseElastic(Vector<IssmDouble>* pUp, Vector<IssmDouble>* pNorth, Vector<IssmDouble>* pEast, Vector<IssmDouble>* pRSLg, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz,int loop,int horiz){/*{{{*/
 
         /*serialized vectors:*/
-        IssmDouble* Sg=NULL;
-
+        IssmDouble* RSLg=NULL;
+        
         IssmDouble  eartharea=0;
         IssmDouble  eartharea_cpu=0;
 
         int         ns,nsmax;
-
+        
         /*Serialize vectors from previous iteration:*/
-        Sg=pSg->ToMPISerial();
+        RSLg=pRSLg->ToMPISerial();
 
         /*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++){
@@ -4446 +4447 @@
         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->SealevelriseGeodetic(pUp,pNorth,pEast,Sg,latitude,longitude,radius,xx,yy,zz,eartharea);
-                }
-                if(i%100==0){
+                        element->SealevelriseGeodetic(pUp,pNorth,pEast,RSLg,latitude,longitude,radius,xx,yy,zz,eartharea,horiz);
+                }
+                if(i%loop==0){
                         pUp->Assemble();
-                        pNorth->Assemble();
-                        pEast->Assemble();
-                }
-        }
-
+                        if (horiz){
+                                pNorth->Assemble();
+                                pEast->Assemble();
+                        }
+                }
+        }
+        
         /*One last time: */
         pUp->Assemble();
-        pNorth->Assemble();
-        pEast->Assemble();
+        if(horiz){
+                pNorth->Assemble();
+                pEast->Assemble();
+        }
+        if(VerboseConvergence())_printf0_("\n");
 
         /*Free ressources:*/
-        xDelete<IssmDouble>(Sg);
-        xDelete<IssmDouble>(latitude);
-        xDelete<IssmDouble>(longitude);
-        xDelete<IssmDouble>(radius);
-        xDelete<IssmDouble>(xx);
-        xDelete<IssmDouble>(yy);
-        xDelete<IssmDouble>(zz);
-}
-/*}}}*/
-IssmDouble FemModel::SealevelriseOceanAverage(Vector<IssmDouble>* Sg) { /*{{{*/
-
-        IssmDouble* Sg_serial=NULL;
+        xDelete<IssmDouble>(RSLg);
+}
+/*}}}*/
+IssmDouble FemModel::SealevelriseOceanAverage(Vector<IssmDouble>* RSLg) { /*{{{*/
+
+        IssmDouble* RSLg_serial=NULL;
         IssmDouble  oceanvalue,oceanvalue_cpu;
         IssmDouble  oceanarea,oceanarea_cpu;
 
         /*Serialize vectors from previous iteration:*/
-        Sg_serial=Sg->ToMPISerial();
+        RSLg_serial=RSLg->ToMPISerial();
 
         /*Initialize:*/
@@ -4488 +4489 @@
                 Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
                 oceanarea_cpu += element->OceanArea();
-                oceanvalue_cpu += element->OceanAverage(Sg_serial);
+                oceanvalue_cpu += element->OceanAverage(RSLg_serial);
         }
         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());
 
         /*Free ressources:*/
-        xDelete<IssmDouble>(Sg_serial);
-
+        xDelete<IssmDouble>(RSLg_serial);
+        
         return oceanvalue/oceanarea;
 }
@@ -4514 +4515 @@
         int*                eplzigzag_counter = NULL;
         int                 eplflip_lock;
-
+        
         HydrologyDCEfficientAnalysis* effanalysis =  new HydrologyDCEfficientAnalysis();
         HydrologyDCInefficientAnalysis* inefanalysis =  new HydrologyDCInefficientAnalysis();
@@ -4521 +4522 @@
         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 +4543 @@
         /*Assemble and serialize*/
         mask->Assemble();
-        serial_mask=mask->ToMPISerial();
-
+        serial_mask=mask->ToMPISerial();        
+        
         xDelete<int>(eplzigzag_counter);
         xDelete<IssmDouble>(serial_rec);
@@ -4585 +4586 @@
         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;
-        Vector<IssmDouble>* active                              = NULL;
-        IssmDouble*         serial_mask = NULL;
-        IssmDouble*         serial_active       = 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 Mask and elementize*/
-        InputUpdateFromVectorx(this,serial_mask,HydrologydcMaskThawedNodeEnum,NodeSIdEnum);
-        xDelete<IssmDouble>(serial_mask);
-        inefanalysis->ElementizeIdsMask(this);
-
-        /*get node mask coherent with element mask*/
-        active=new Vector<IssmDouble>(nodes->NumberOfNodes(HydrologyDCInefficientAnalysisEnum));
-        for (int i=0;i<elements->Size();i++){
-                Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
-                inefanalysis->HydrologyIdsGetActive(active,element);
-        }
-
-        /*Assemble and serialize*/
-        active->Assemble();
-        serial_active=active->ToMPISerial();
-        delete active;
-
-        /*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_active[node->Sid()]==1.){
-                                node->Activate();
-                                if(!node->IsClone()) counter++;
-                        }
-                        else{
-                                node->Deactivate();
-                        }
-                }
-        }
-
-        xDelete<IssmDouble>(serial_active);
-        delete inefanalysis;
-        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*/
@@ -4706 +4632 @@
         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;
@@ -4791 +4717 @@
 }
 /*}}}*/
-#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: */
@@ -4806 +4732 @@
         /*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));
@@ -4826 +4752 @@
         size_t* poutputbuffersize;
 
-
+        
         /*Before we delete the profiler, report statistics for this run: */
         profiler->Stop(TOTAL);  //final tagging
@@ -4839 +4765 @@
                                 );
         _printf0_("\n");
-
+        
         /*Before we close the output file, recover the buffer and size:*/
         outputbufferparam = xDynamicCast<GenericParam<char**>*>(this->parameters->FindParamObject(OutputBufferPointerEnum));
@@ -4879 +4805 @@
 
         /*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);
@@ -4887 +4813 @@
 }/*}}}*/
 #endif
+
 
 #if defined(_HAVE_BAMG_) && !defined(_HAVE_ADOLC_)
@@ -4920 +4847 @@
                 /*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);
@@ -4927 +4854 @@
                 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);
@@ -4947 +4874 @@
                 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*/
@@ -4983 +4910 @@
         /*Create bamg data structures for bamg*/
         this->amrbamg = new AmrBamg();
-
+        
         /*Get amr parameters*/
         this->parameters->FindParam(&hmin,AmrHminEnum);
@@ -5007 +4934 @@
 
         /*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);
         }
@@ -5021 +4948 @@
 
         if(!hmaxvertices) _error_("hmaxvertices is NULL!\n");
-
+        
         /*Intermediaries*/
         int numberofvertices                     = this->vertices->NumberOfVertices();
@@ -5028 +4955 @@
 
         switch(levelset_type){
-                case MaskGroundediceLevelsetEnum:
+                case MaskGroundediceLevelsetEnum: 
                         threshold       = this->amrbamg->groundingline_distance;
                         resolution      = this->amrbamg->groundingline_resolution;
@@ -5042 +4969 @@
         this->GetVerticeDistanceToZeroLevelSet(&verticedistance,levelset_type);
         if(!verticedistance) _error_("verticedistance is NULL!\n");
-
+        
         /*Fill hmaxVertices*/
         for(int i=0;i<numberofvertices;i++){
@@ -5056 +4983 @@
 /*}}}*/
 void FemModel::GethmaxVerticesFromEstimators(IssmDouble* hmaxvertices,int errorestimator_type){/*{{{*/
-
+   
         if(!hmaxvertices) _error_("hmaxvertices is NULL!\n");
 
@@ -5064 +4991 @@
         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;
@@ -5077 +5004 @@
         /*Fill variables*/
         switch(errorestimator_type){
-                case ThicknessErrorEstimatorEnum:
+                case ThicknessErrorEstimatorEnum: 
                         threshold               = this->amrbamg->thicknesserror_threshold;
                         groupthreshold  = this->amrbamg->thicknesserror_groupthreshold;
@@ -5102 +5029 @@
         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++){
@@ -5121 +5048 @@
                 error_vertices[v2]+=error_elements[i];
                 error_vertices[v3]+=error_elements[i];
-        }
+        }       
 
         /*Fill hmaxvertices with the criteria*/
@@ -5135 +5062 @@
                         }
                 }
-                /*Now, fill the hmaxvertices if requested*/
+                /*Now, fill the hmaxvertices if requested*/       
                 if(refine){
                         for(int j=0;j<elementswidth;j++){
@@ -5165 +5092 @@
         /*Output*/
         IssmDouble* verticedistance;
-
+        
         /*Intermediaries*/
    int numberofvertices       = this->vertices->NumberOfVertices();
@@ -5177 +5104 @@
         /*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);
 
@@ -5187 +5114 @@
                 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*/
@@ -5223 +5150 @@
         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.");
@@ -5242 +5169 @@
                 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;
@@ -5263 +5190 @@
 /*}}}*/
 void FemModel::InitializeAdaptiveRefinementNeopz(void){/*{{{*/
-
+        
         /*Define variables*/
         int my_rank                                                                             = IssmComm::GetRank();
@@ -5272 +5199 @@
         IssmDouble* z                                                                   = NULL;
         int* elements                                                                   = NULL;
-
+        
         /*Initialize field as NULL for now*/
         this->amr = NULL;
@@ -5280 +5207 @@
         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);
@@ -5298 +5225 @@
         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);
         }
@@ -5319 +5246 @@
         /*Output*/
         IssmDouble* elementdistance;
-
+        
         /*Intermediaries*/
    int numberofelements                                                 = this->elements->NumberOfElements();
@@ -5328 +5255 @@
         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);
 
@@ -5335 +5262 @@
       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];
@@ -5345 +5272 @@
                 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*/