Changeset 21574

Timestamp:

02/22/17 11:34:34 (8 years ago)

Author:

tsantos

Message:

BUG: fixed a bug in CreateConstaints (adaptive mesh refinement method): bad use of Vector<IssmDouble>. NOTE: AMR not working yet

File:

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

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

@@ -2954 +2954 @@
 /*}}}*/
 void FemModel::ReMesh(void){/*{{{*/
-        
+
         /*Variables*/
         IssmDouble *newx                        = NULL;
@@ -2973 +2973 @@
         this->ElementsAndVerticesPartitioning(newnumberofvertices,newnumberofelements,elementswidth,newelementslist,&my_elements,&my_vertices);
 
-        if(this->loads->Size()!=0){
-                _error_("not supported yet");
-        }
+        if(this->loads->Size()!=0) _error_("not supported yet");
 
         /*Create vertices*/
@@ -3009 +3007 @@
                 this->UpdateElements(newnumberofelements,newelementslist,my_elements,nodecounter,i,new_elements);
 
-
-                if(IssmComm::GetRank()==2){
-                        _printf_("New constraints: CPU 1\n");
-                        //this->constraints->DeepEcho();
-                        new_constraints->DeepEcho();
-                        printarray(my_vertices,newnumberofvertices,1);
-                        _error_("stop!");
-                }
-
-
                 if(new_nodes->Size()) nodecounter=new_nodes->MaximumId();
                 constraintcounter = new_constraints->NumberOfConstraints();
@@ -3069 +3057 @@
         /*Finally: interpolate all inputs and insert them into the new elements.*/
         this->InterpolateInputs(new_vertices,new_elements);
-
-        if(IssmComm::GetRank()==2){
-                _printf_("New constraints: CPU 1\n");
-                //this->constraints->DeepEcho();
-                new_constraints->DeepEcho();
-                _error_("stop!");
-        }
-
 
         /*Delete old structure and set new pointers*/
@@ -3631 +3611 @@
         Vector<IssmDouble>* vspcvyflag  = new Vector<IssmDouble>(numberofnodes_analysistype);
         
-        for(int i=0;i<numberofnodes_analysistype;i++){ 
-                vspcvx->SetValue(i,0.,INS_VAL);
-                vspcvy->SetValue(i,0.,INS_VAL);
-                vspcvxflag->SetValue(i,0.,INS_VAL);
-                vspcvyflag->SetValue(i,0.,INS_VAL);
-        }
-
         for(int i=0;i<this->constraints->Size();i++){
                 SpcStatic* spc                  = xDynamicCast<SpcStatic*>(this->constraints->GetObjectByOffset(i));
@@ -3645 +3618 @@
                 int nodeindex                   = node-1;
                 
-                if(IssmComm::GetRank()==0)      
-                  {//itapopo
-                 if(dof==0) {//vx
-                        if(IssmComm::GetRank()==0) _printf_("Node: " << nodeindex << "\tdof: " << dof << "\tspcvalue: " << spcvalue << "\n") ;
+                if(dof==0) {//vx
                         vspcvx->SetValue(nodeindex,spcvalue,INS_VAL);
                         vspcvxflag->SetValue(nodeindex,1,INS_VAL);
@@ -3656 +3626 @@
                         vspcvyflag->SetValue(nodeindex,1,INS_VAL);
                 }
-
-                  }//itapopo
-        }
-
-        #ifdef _HAVE_PETSC_
-                _printf_("Tem PETSC\n");
-        #endif
+        }
 
         /*Assemble*/
@@ -3670 +3634 @@
         vspcvyflag->Assemble();
 
-        IssmDouble pvalue1,pvalue2;
-        vspcvx->GetValue(&pvalue1,0);
-        vspcvxflag->GetValue(&pvalue2,0);
-        if(IssmComm::GetRank()==0) _printf_("Node: " << 0 << "\tspcvalue: " << pvalue1 << "\tflag: " << pvalue2 << "\n") ;
-        vspcvx->GetValue(&pvalue1,5);
-        vspcvxflag->GetValue(&pvalue2,5);
-        if(IssmComm::GetRank()==0) _printf_("Node: " << 5 << "\tspcvalue: " << pvalue1 << "\tflag: " << pvalue2 << "\n") ;
-        
-        
-        
         /*Serialize*/
         spcvx            = vspcvx->ToMPISerial();
@@ -3702 +3656 @@
         InterpFromMeshToMesh2dx(&newspcvyflag,elementslist,x,y,nods_data,nels_data,spcvyflag,M_data,N_data,newx,newy,N_interp,NULL);
 
-        if(IssmComm::GetRank()==1){
-                _printf_("CPU 1:\n");
-                _printf_("Old spcvx:\n");
-                printarray(spcvx,nods_data,1);
-                _printf_("Old spcvxflag:\n");
-                printarray(spcvxflag,nods_data,1);
-                _error_("Stop!!");
-        }
-
-
-
-
-
         int count                                       = 0;
         IssmDouble eps                          = 1.e-8;
@@ -3720 +3661 @@
         /*Now, insert the interpolated constraints in the data set (constraints)*/
         for(int i=0;i<newnumberofvertices;i++){
-                if(my_vertices[i])
-                /*spcvx*/
-                if(!xIsNan<IssmDouble>(newspcvx[i]) && newspcvxflag[i]>(1-eps)){
-                        constraints->AddObject(new SpcStatic(constraintcounter+count+1,nodecounter+i+1,0,newspcvx[i],StressbalanceAnalysisEnum));
-                        //add count'th spc, on node i+1, setting dof 1 to vx.
-                        count++;
+                if(my_vertices[i]){ 
+                        /*spcvx*/
+                        if(!xIsNan<IssmDouble>(newspcvx[i]) && newspcvxflag[i]>(1-eps)){
+                                constraints->AddObject(new SpcStatic(constraintcounter+count+1,nodecounter+i+1,0,newspcvx[i],StressbalanceAnalysisEnum));
+                                //add count'th spc, on node i+1, setting dof 1 to vx.
+                                count++;
+                        }
                 }
         }
         count=0;
         for(int i=0;i<newnumberofvertices;i++){
-                if(my_vertices[i])
-                /*spcvy*/
-                if(!xIsNan<IssmDouble>(newspcvy[i]) && newspcvyflag[i]>(1-eps) ){
-                        constraints->AddObject(new SpcStatic(constraintcounter+count+1,nodecounter+i+1,1,newspcvy[i],StressbalanceAnalysisEnum)); 
-                        //add count'th spc, on node i+1, setting dof 1 to vx.
-                        count++;
+                if(my_vertices[i]){
+                        /*spcvy*/
+                        if(!xIsNan<IssmDouble>(newspcvy[i]) && newspcvyflag[i]>(1-eps) ){
+                                constraints->AddObject(new SpcStatic(constraintcounter+count+1,nodecounter+i+1,1,newspcvy[i],StressbalanceAnalysisEnum)); 
+                                //add count'th spc, on node i+1, setting dof 1 to vx.
+                                count++;
+                        }
                 }
         }