Changeset 23046

Timestamp:

08/02/18 16:44:35 (7 years ago)

Author:

jdquinn

Message:

CHG: Reverted mistaken commit

Location:

issm/trunk-jpl

Files:

• jenkins/macosx_pine-island (modified) (2 diffs)
• jenkins/macosx_pine-island_static (modified) (3 diffs)
• src/c/classes/FemModel.cpp (modified) (139 diffs)
• src/c/modules/InterpFromGridToMeshx/InterpFromGridToMeshx.cpp (modified) (6 diffs)
• src/wrappers/InterpFromGridToMesh/InterpFromGridToMesh.js (modified) (16 diffs)

issm/trunk-jpl/jenkins/macosx_pine-island

@@ -5 +5 @@
 
 #MATLAB path
-MATLAB_PATH="/Applications/MATLAB_R2017b.app/"
+MATLAB_PATH="/Applications/MATLAB_R2015b.app/"
 
-#ISSM CONFIGURATION 
+#ISSM CONFIGURATION
 ISSM_CONFIG='--prefix=$ISSM_DIR \
         --with-matlab-dir=$MATLAB_PATH \
@@ -54 +54 @@
 #by Matlab and runme.m
 #ex: "'id',[101 102 103]"
-##                           FS                     
+##                           FS
 PYTHON_NROPTIONS=""
 MATLAB_NROPTIONS="'exclude',[701,702,703,435,IdFromString('Dakota')]"

issm/trunk-jpl/jenkins/macosx_pine-island_static

@@ -5 +5 @@
 
 #MATLAB path
-MATLAB_PATH="/Applications/MATLAB_R2017b.app/"
+MATLAB_PATH="/Applications/MATLAB_R2015b.app/"
 
-#ISSM CONFIGURATION 
+#ISSM CONFIGURATION
 ISSM_CONFIG='--prefix=$ISSM_DIR \
         --disable-static \
@@ -23 +23 @@
         --with-m1qn3-dir=$ISSM_DIR/externalpackages/m1qn3/install \
         --with-math77-dir=$ISSM_DIR/externalpackages/math77/install \
-        --with-fortran-lib="/usr/local/gfortran/lib/libgfortran.a /usr/local/gfortran/lib/libquadmath.a /usr/local/gfortran/lib/gcc/x86_64-apple-darwin15/6.1.0/libgcc.a" \
+        --with-fortran-lib="/usr/local/gfortran/lib/libgfortran.a /usr/local/gfortran/lib/libquadmath.a /usr/local/gfortran/lib/gcc/x86_64-apple-darwin14/5.2.0/libgcc.a" \
         --with-numthreads=4'
 
@@ -61 +61 @@
 #by Matlab and runme.m
 #ex: "'id',[101 102 103]"
-##                           bamg mesh   FS                     
+##                           bamg mesh   FS
 #PYTHON_NROPTIONS=""
 #MATLAB_NROPTIONS="'exclude',[119,243,514,701,702,703,435,IdFromString('Dakota')]"

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

@@ -1464 +1464 @@
         /*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;
 
@@ -1488 +1488 @@
         /*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;
 
@@ -1512 +1512 @@
         /*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;
 
@@ -1536 +1536 @@
         /*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;
 
@@ -1560 +1560 @@
         /*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;
 
@@ -1584 +1584 @@
         /*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;
 
@@ -1629 +1629 @@
                         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;
@@ -1965 +1965 @@
                                                                 IssmDouble* values    = xNewZeroInit<IssmDouble>(nlines*ncols);
                                                                 IssmDouble* allvalues = xNew<IssmDouble>(nlines*ncols);
-                                                                
+
                                                                 /*Fill-in matrix*/
                                                                 for(int j=0;j<elements->Size();j++){
@@ -1974 +1974 @@
                                                                 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);
@@ -2075 +2075 @@
                 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;
         }
 
@@ -2212 +2212 @@
                         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;
                 }
@@ -2403 +2403 @@
         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: */
@@ -2416 +2416 @@
         IssmDouble         *surface = NULL;
         IssmDouble         *bed     = NULL;
-                        
+
         if(VerboseSolution()) _printf0_("   updating vertices positions\n");
 
@@ -2455 +2455 @@
 
 /*AMR*/
-#if !defined(_HAVE_ADOLC_) 
+#if !defined(_HAVE_ADOLC_)
 void FemModel::ReMesh(void){/*{{{*/
 
@@ -2465 +2465 @@
         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
@@ -2471 +2471 @@
         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){
@@ -2484 +2484 @@
                 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);
@@ -2515 +2515 @@
 
                 /*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);
@@ -2545 +2545 @@
                 //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);
@@ -2562 +2562 @@
 
                 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
                 }
@@ -2606 +2606 @@
         /*Insert bedrock from mismip+ setup*/
         /*This was used to Misomip project/simulations*/
-        
+
         if(VerboseSolution())_printf0_("        call Mismip bedrock adjust module\n");
 
@@ -2623 +2623 @@
                         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);
@@ -2636 +2636 @@
 
         if(VerboseSolution())_printf0_("        call adjust base and thickness module\n");
-        
+
         int     numvertices = this->GetElementsWidth();
    IssmDouble rho_water,rho_ice,density,base_float;
@@ -2648 +2648 @@
         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);
@@ -2660 +2660 @@
                         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;
                 }
 
@@ -2677 +2677 @@
                 element->AddInput(BaseEnum,&b[0],P1Enum);
         }
-        
+
    /*Delete*/
    xDelete<IssmDouble>(phi);
@@ -2705 +2705 @@
         Vector<IssmDouble>* input_interpolations        = NULL;
         IssmDouble* input_interpolations_serial = NULL;
-   int* pos                                                                                             = NULL; 
+   int* pos                                                                                             = NULL;
         IssmDouble value                                                                        = 0;
 
@@ -2725 +2725 @@
                         P0input_interp = xNew<int>(numP0inputs);
                         P1input_enums  = xNew<int>(numP1inputs);
-                        P1input_interp = xNew<int>(numP1inputs);        
+                        P1input_interp = xNew<int>(numP1inputs);
                 }
                 numP0inputs = 0;
@@ -2757 +2757 @@
                 }
         }
-        
-        /*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);
@@ -2764 +2764 @@
         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++){
@@ -2780 +2780 @@
                         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);
                 }
         }
@@ -2790 +2790 @@
         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*/
@@ -2810 +2810 @@
 /*}}}*/
 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
@@ -2826 +2826 @@
         int* P1input_enums                      = NULL;
         int* P1input_interp                     = NULL;
-        IssmDouble* values                      = NULL; 
+        IssmDouble* values                      = NULL;
    IssmDouble* vector           = NULL;
         IssmDouble* x                                   = NULL;//global, entire old mesh
@@ -2838 +2838 @@
         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)*/
@@ -2871 +2871 @@
 
         /*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;
@@ -2899 +2899 @@
                 }
         }
-        
+
         /*Cleanup*/
         xDelete<IssmDouble>(P0inputs);
@@ -2938 +2938 @@
 
         if(!this->elements || !this->vertices || !this->results || !this->parameters) return;
-         
+
         parameters->FindParam(&step,StepEnum);
         parameters->FindParam(&time,TimeEnum);
@@ -2956 +2956 @@
         this->results->AddResult(new GenericExternalResult<IssmDouble*>(this->results->Size()+1,MeshYEnum,
                                         y,numberofvertices,1,step,time));
-        
+
         /*Cleanup*/
         xDelete<IssmDouble>(x);
@@ -3017 +3017 @@
    /*Assemble*/
         vmasklevelset->Assemble();
-        
+
         /*Serialize and set output*/
         (*pmasklevelset)=vmasklevelset->ToMPISerial();
@@ -3031 +3031 @@
 
         /*newelementslist is in Matlab indexing*/
-        
+
         /*Creating connectivity table*/
         int* connectivity=NULL;
@@ -3042 +3042 @@
                         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;
@@ -3058 +3058 @@
                         newvertex->connectivity=connectivity[i];
                         newvertex->clone=false;//itapopo check this
-                        vertices->AddObject(newvertex); 
-                } 
+                        vertices->AddObject(newvertex);
+                }
         }
 
@@ -3086 +3086 @@
                         }
                         else newtria->element_type_list=NULL;
-                        
+
                         /*Element hook*/
                         int matpar_id=newnumberofelements+1; //retrieve material parameter id (last pointer in femodel->materials)
@@ -3092 +3092 @@
                         /*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;
@@ -3104 +3104 @@
                         /*Clean up*/
                         xDelete<int>(vertex_ids);
-                        elements->AddObject(newtria);   
-                } 
+                        elements->AddObject(newtria);
+                }
         }
 
@@ -3115 +3115 @@
         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
 }
 /*}}}*/
@@ -3129 +3129 @@
         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;
@@ -3138 +3138 @@
                         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;
@@ -3159 +3159 @@
                         /*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]);
@@ -3174 +3174 @@
         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();
@@ -3180 +3180 @@
         IssmDouble* x                   = NULL;
         IssmDouble* y                   = NULL;
-        IssmDouble* z                   = NULL; 
+        IssmDouble* z                   = NULL;
         int* elementslist       = NULL;
         int* elem_vertices      = NULL;
@@ -3189 +3189 @@
         /*Get vertices coordinates*/
         VertexCoordinatesx(&x,&y,&z,femmodel_vertices,false) ;
-        
+
         /*Get element vertices*/
         elem_vertices                           = xNew<int>(elementswidth);
@@ -3204 +3204 @@
         vid3->SetValue(element->sid,elem_vertices[2],INS_VAL);
    }
-                
+
         /*Assemble*/
    vid1->Assemble();
@@ -3214 +3214 @@
    id2 = vid2->ToMPISerial();
         id3 = vid3->ToMPISerial();
-        
+
         /*Construct elements list*/
         elementslist=xNew<int>(numberofelements*elementswidth);
@@ -3223 +3223 @@
                 elementslist[elementswidth*i+2] = reCast<int>(id3[i])+1; //InterpMesh wants Matlab indexing
         }
-        
+
         /*Assign pointers*/
         *px                             = x;
@@ -3244 +3244 @@
         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)
@@ -3251 +3251 @@
         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
@@ -3261 +3261 @@
         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();
@@ -3277 +3277 @@
         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));
@@ -3284 +3284 @@
                 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;
@@ -3291 +3291 @@
         *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*/
@@ -3302 +3302 @@
         /*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
@@ -3328 +3328 @@
         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();
@@ -3336 +3336 @@
         /*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");
@@ -3343 +3343 @@
                 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
@@ -3366 +3366 @@
                                                                 spc,numberofvertices,numberofcols,
                                                                 newx,newy,newnumberofvertices,NULL);
-        
+
         /*Now, insert the interpolated constraints in the data set (constraints)*/
         count=0;
@@ -3383 +3383 @@
                 /*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++;
@@ -3440 +3440 @@
         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++){
@@ -3466 +3466 @@
                 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);
@@ -3476 +3476 @@
         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
                         }
@@ -3495 +3495 @@
         /*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);
@@ -3516 +3516 @@
    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++){
@@ -3564 +3564 @@
         /*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];
@@ -3589 +3589 @@
 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*/
 
@@ -3609 +3609 @@
         /*Get smoothed deviatoric stress tensor*/
         this->SmoothedDeviatoricStressTensor(&smoothedtauxx,&smoothedtauyy,&smoothedtauxy);
-        
+
         /*Integrate the error over elements*/
    for(int i=0;i<this->elements->Size();i++){
@@ -3617 +3617 @@
       element->GetInputListOnVertices(tauxy,DeviatoricStressxyEnum);
       element->GetVerticesSidList(elem_vertices);
-                
+
                 /*Integrate*/
                 element->GetVerticesCoordinates(&xyz_list);
@@ -3632 +3632 @@
                                 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)
@@ -3648 +3648 @@
    /*Serialize and set output*/
    (*pelementerror)=velementerror->ToMPISerial();
-        
+
         /*Cleanup*/
         xDelete<IssmDouble>(smoothedtauxx);
@@ -3692 +3692 @@
       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);
@@ -3760 +3760 @@
    /*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++){
@@ -3846 +3846 @@
         IssmDouble* x                                   = NULL;
         IssmDouble* y                                   = NULL;
-        IssmDouble* z                                   = NULL; 
+        IssmDouble* z                                   = NULL;
         IssmDouble* xyz_list                    = NULL;
         IssmDouble x1,y1,x2,y2,x3,y3;
@@ -3855 +3855 @@
       //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];
@@ -3888 +3888 @@
                 _error_("level set type not implemented yet!");
         }
-        
+
         /*Outputs*/
         IssmDouble* zerolevelset_points                 = NULL;
         int npoints                                                                             = 0;
-        
+
         /*Intermediaries*/
         int elementswidth                       = this->GetElementsWidth();
@@ -3905 +3905 @@
         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++){
@@ -3912 +3912 @@
                 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.;
@@ -3950 +3950 @@
                 }
         }
-        
+
         /*Assign outputs*/
         numberofpoints                          = npoints;
@@ -3990 +3990 @@
         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: */
 
@@ -4083 +4083 @@
 
         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());
@@ -4105 +4105 @@
                 }
         }
-        
+
         /*One last time: */
         pUp->Assemble();
@@ -4124 +4124 @@
 
         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));
@@ -4152 +4152 @@
                 }
         }
-        
+
         /*One last time: */
         pUp->Assemble();
@@ -4214 +4214 @@
         }
         if(VerboseConvergence())_printf0_("\n");
-                
+
         /*One last time: */
         pRSLgi->Assemble();
@@ -4232 +4232 @@
         /*serialized vectors:*/
         IssmDouble* RSLg_old=NULL;
-        
+
         IssmDouble  eartharea=0;
         IssmDouble  eartharea_cpu=0;
 
         int         ns,nsmax;
-        
+
         /*Serialize vectors from previous iteration:*/
         RSLg_old=pRSLg_old->ToMPISerial();
@@ -4249 +4249 @@
                 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());
@@ -4267 +4267 @@
         }
         if(verboseconvolution)if(VerboseConvergence())_printf0_("\n");
-        
+
         /*Free ressources:*/
         xDelete<IssmDouble>(RSLg_old);
@@ -4355 +4355 @@
         /*Free ressources:*/
         xDelete<IssmDouble>(RSLg_old);
-        
-        
+
+
 }
 /*}}}*/
@@ -4363 +4363 @@
         /*serialized vectors:*/
         IssmDouble* RSLg=NULL;
-        
+
         IssmDouble  eartharea=0;
         IssmDouble  eartharea_cpu=0;
 
         int         ns,nsmax;
-        
+
         /*Serialize vectors from previous iteration:*/
         RSLg=pRSLg->ToMPISerial();
@@ -4374 +4374 @@
         /*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++){
@@ -4402 +4402 @@
                 }
         }
-        
+
         /*One last time: */
         pUp->Assemble();
@@ -4436 +4436 @@
         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());
@@ -4442 +4442 @@
         /*Free ressources:*/
         xDelete<IssmDouble>(RSLg_serial);
-        
+
         return oceanvalue/oceanarea;
 }
@@ -4458 +4458 @@
         int*                eplzigzag_counter = NULL;
         int                 eplflip_lock;
-        
+
         HydrologyDCEfficientAnalysis* effanalysis =  new HydrologyDCEfficientAnalysis();
         HydrologyDCInefficientAnalysis* inefanalysis =  new HydrologyDCInefficientAnalysis();
@@ -4465 +4465 @@
         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));
@@ -4486 +4486 @@
         /*Assemble and serialize*/
         mask->Assemble();
-        serial_mask=mask->ToMPISerial();        
-        
+        serial_mask=mask->ToMPISerial();
+
         xDelete<int>(eplzigzag_counter);
         xDelete<IssmDouble>(serial_rec);
@@ -4529 +4529 @@
         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;
@@ -4650 +4650 @@
         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;
@@ -4735 +4735 @@
 }
 /*}}}*/
-#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: */
@@ -4750 +4750 @@
         /*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));
@@ -4770 +4770 @@
         size_t* poutputbuffersize;
 
-        
+
         /*Before we delete the profiler, report statistics for this run: */
         profiler->Stop(TOTAL);  //final tagging
@@ -4783 +4783 @@
                                 );
         _printf0_("\n");
-        
+
         /*Before we close the output file, recover the buffer and size:*/
         outputbufferparam = xDynamicCast<GenericParam<char**>*>(this->parameters->FindParamObject(OutputBufferPointerEnum));
@@ -4823 +4823 @@
 
         /*Open output file once for all and add output file descriptor to parameters*/
-        output_fid=open_memstream(&outputbuffer,&outputsize); 
-        if(output_fid==NULL||output_fid==0)_error_("could not initialize output stream");
+        output_fid=open_memstream(&outputbuffer,&outputsize);
+        if(output_fid==NULL)_error_("could not initialize output stream");
         this->parameters->SetParam(output_fid,OutputFilePointerEnum);
         this->parameters->AddObject(new GenericParam<char**>(&outputbuffer,OutputBufferPointerEnum));
@@ -4865 +4865 @@
                 /*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);
@@ -4872 +4872 @@
                 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);
@@ -4892 +4892 @@
                 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*/
@@ -4928 +4928 @@
         /*Create bamg data structures for bamg*/
         this->amrbamg = new AmrBamg();
-        
+
         /*Get amr parameters*/
         this->parameters->FindParam(&hmin,AmrHminEnum);
@@ -4952 +4952 @@
 
         /*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);
         }
@@ -4966 +4966 @@
 
         if(!hmaxvertices) _error_("hmaxvertices is NULL!\n");
-        
+
         /*Intermediaries*/
         int numberofvertices                     = this->vertices->NumberOfVertices();
@@ -4973 +4973 @@
 
         switch(levelset_type){
-                case MaskGroundediceLevelsetEnum: 
+                case MaskGroundediceLevelsetEnum:
                         threshold       = this->amrbamg->groundingline_distance;
                         resolution      = this->amrbamg->groundingline_resolution;
@@ -4987 +4987 @@
         this->GetVerticeDistanceToZeroLevelSet(&verticedistance,levelset_type);
         if(!verticedistance) _error_("verticedistance is NULL!\n");
-        
+
         /*Fill hmaxVertices*/
         for(int i=0;i<numberofvertices;i++){
@@ -5001 +5001 @@
 /*}}}*/
 void FemModel::GethmaxVerticesFromEstimators(IssmDouble* hmaxvertices,int errorestimator_type){/*{{{*/
-   
+
         if(!hmaxvertices) _error_("hmaxvertices is NULL!\n");
 
@@ -5009 +5009 @@
         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;
@@ -5022 +5022 @@
         /*Fill variables*/
         switch(errorestimator_type){
-                case ThicknessErrorEstimatorEnum: 
+                case ThicknessErrorEstimatorEnum:
                         threshold               = this->amrbamg->thicknesserror_threshold;
                         groupthreshold  = this->amrbamg->thicknesserror_groupthreshold;
@@ -5047 +5047 @@
         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++){
@@ -5066 +5066 @@
                 error_vertices[v2]+=error_elements[i];
                 error_vertices[v3]+=error_elements[i];
-        }       
+        }
 
         /*Fill hmaxvertices with the criteria*/
@@ -5080 +5080 @@
                         }
                 }
-                /*Now, fill the hmaxvertices if requested*/       
+                /*Now, fill the hmaxvertices if requested*/
                 if(refine){
                         for(int j=0;j<elementswidth;j++){
@@ -5110 +5110 @@
         /*Output*/
         IssmDouble* verticedistance;
-        
+
         /*Intermediaries*/
    int numberofvertices       = this->vertices->NumberOfVertices();
@@ -5122 +5122 @@
         /*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);
 
@@ -5132 +5132 @@
                 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*/
@@ -5168 +5168 @@
         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.");
@@ -5187 +5187 @@
                 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;
@@ -5208 +5208 @@
 /*}}}*/
 void FemModel::InitializeAdaptiveRefinementNeopz(void){/*{{{*/
-        
+
         /*Define variables*/
         int my_rank                                                                             = IssmComm::GetRank();
@@ -5217 +5217 @@
         IssmDouble* z                                                                   = NULL;
         int* elements                                                                   = NULL;
-        
+
         /*Initialize field as NULL for now*/
         this->amr = NULL;
@@ -5225 +5225 @@
         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);
@@ -5243 +5243 @@
         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);
         }
@@ -5264 +5264 @@
         /*Output*/
         IssmDouble* elementdistance;
-        
+
         /*Intermediaries*/
    int numberofelements                                                 = this->elements->NumberOfElements();
@@ -5273 +5273 @@
         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);
 
@@ -5280 +5280 @@
       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];
@@ -5290 +5290 @@
                 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/modules/InterpFromGridToMeshx/InterpFromGridToMeshx.cpp

@@ -1 +1 @@
 /*!\file:  InterpFromGridToMeshx.cpp
  * \brief  "c" core code for interpolating values from a structured grid.
- */ 
+ */
 
 #ifdef HAVE_CONFIG_H
@@ -24 +24 @@
         double* y=NULL;
         int     i;
-
-        // TEST
-        _printf_("\r      N: "<<N<<" M:"<<M<<"\n");
 
         /*Some checks on arguments: */
@@ -151 +148 @@
                          *    |         |     |
                          * y1 x---------+-----x Q21
-                         *    x1                 x2       
+                         *    x1                 x2
                          *
                          */
@@ -225 +222 @@
         /*split the rectangle in 2 triangle and
          * use Lagrange P1 interpolation
-         *       
+         *
          *   +3----+2,3' Q12----Q22
          *   |    /|     |    /|
@@ -282 +279 @@
         _assert_(x<=x2 && x>=x1 && y<=y2 && y>=y1);
 
-        return 
+        return
           +Q11*(x2-x)*(y2-y)/((x2-x1)*(y2-y1))
           +Q21*(x-x1)*(y2-y)/((x2-x1)*(y2-y1))
@@ -301 +298 @@
          *    |        |         |
          * y1 x--------x---------x Q21
-         *    x1       xm        x2 
+         *    x1       xm        x2
          *
          */

issm/trunk-jpl/src/wrappers/InterpFromGridToMesh/InterpFromGridToMesh.js

@@ -6 +6 @@
  * Usage:
  *      var data_mesh=InterpFromGridToMesh(xIn,yIn,dataIn,xMeshIn,yMeshIn,defaultValue,interpType);\
- * 
+ *
  *      xIn,yIn                                         : coordinates of matrix data. (x and y must be in increasing order)
  *      dataIn                                          : matrix holding the data to be interpolated onto the mesh
@@ -51 +51 @@
         var yMesh                       = {};
         //}}}
-        
+
 
         /*
@@ -57 +57 @@
         */
         //{{{
-        
+
         /*
                 Input
@@ -68 +68 @@
         dxHeap.set(new Uint8Array(dx.buffer));
         x                       = dxHeap.byteOffset;
-        
+
         dy                      = new Float64Array(yIn);
         ny                      = dy.length * dy.BYTES_PER_ELEMENT;
@@ -75 +75 @@
         dyHeap.set(new Uint8Array(dy.buffer));
         y                       = dyHeap.byteOffset;
-        
+
         ddata           = new Float64Array(dataIn);
         ndata           = ddata.length * ddata.BYTES_PER_ELEMENT;
@@ -82 +82 @@
         ddataHeap.set(new Uint8Array(ddata.buffer));
         data            = ddataHeap.byteOffset;
-        
+
         dxMesh          = new Float64Array(xMeshIn);
         nxMesh          = dxMesh.length * dxMesh.BYTES_PER_ELEMENT;
@@ -89 +89 @@
         dxMeshHeap.set(new Uint8Array(dxMesh.buffer));
         xMesh           = dxMeshHeap.byteOffset;
-        
+
         dyMesh          = new Float64Array(yMeshIn);
         nyMesh          = dyMesh.length * dyMesh.BYTES_PER_ELEMENT;
@@ -96 +96 @@
         dyMeshHeap.set(new Uint8Array(dyMesh.buffer));
         yMesh           = dyMeshHeap.byteOffset;
-        
+
         nods            = xMeshIn.length;
         meshNumRows     = xMeshIn.length;
-        console.log('InterpFromGridToMesh.js: meshNumRows => ' + meshNumRows);
-        
-        
+
+
         /*
                 Retrieve interpolation type
@@ -112 +111 @@
         }
         //}}}
-        
+
         /*
                 Output
@@ -118 +117 @@
         pdataMesh = Module._malloc(4);
         //}}}
-        
+
         //}}}
-        
-        
+
+
         /*
                 Declare InterpFromGridToMesh module
@@ -127 +126 @@
         //{{{
         InterpFromGridToMeshModule = Module.cwrap(
-                'InterpFromGridToMeshModule', 
-                'number', 
+                'InterpFromGridToMeshModule',
+                'number',
                 [
                         'number', // output : pdataMesh
@@ -145 +144 @@
         );
         //}}}
-        
-        
+
+
         /*
                 Call InterpFromGridToMesh module
@@ -152 +151 @@
         //{{{
         InterpFromGridToMeshModule(
-                pdataMesh, 
-                x, 
-                y, 
-                data, 
-                xMesh, 
-                yMesh, 
-                defaultValue, 
+                pdataMesh,
+                x,
+                y,
+                data,
+                xMesh,
+                yMesh,
+                defaultValue,
                 nods,
                 dataNumRowsIn,
@@ -166 +165 @@
         );
         //}}}
-        
-        
+
+
         /*
                 Dynamic copying from heap
@@ -175 +174 @@
         dataMesh        = Module.HEAPF64.slice(dataMeshPtr / 8, dataMeshPtr / 8 + nods);
         //}}}
-        
-        
+
+
         /*
                 Free resources
@@ -183 +182 @@
         Module._free(pdataMesh);
         //}}}
-        
-        
+
+
         return dataMesh;
 }