Changeset 15252 for issm/trunk-jpl/src/c/classes/Elements/Tria.cpp

Timestamp:

06/13/13 09:45:43 (12 years ago)

Author:

Mathieu Morlighem

Message:

CHG: moving mask transformation from UpdateFromSolution to GetMask

File:

• issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (modified) (6 diffs)

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

@@ -6407 +6407 @@
 
         /*Intermediaries*/
-        const int   numdof         = NDOF1 *NUMVERTICES;
-        int        *doflist        = NULL;
+        const int   numdof   = NDOF1 *NUMVERTICES;
+        int        *doflist  = NULL;
         bool        converged;
         bool        isefficientlayer;
-        IssmDouble  activeEpl[numdof],activate[numdof]={0.0};
         IssmDouble  values[numdof];
         IssmDouble  residual[numdof];
@@ -6433 +6432 @@
 
         /*Get inputs*/
-        if(isefficientlayer){
-                GetInputListOnVertices(&activeEpl[0],HydrologydcMaskEplactiveEnum);
-        }
         if(converged){
                 this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
                 this->parameters->FindParam(&kmax,HydrologySedimentKmaxEnum);
                 this->parameters->FindParam(&penalty_factor,HydrologydcPenaltyFactorEnum);
-                
                 
                 kappa=kmax*pow(10.,penalty_factor);
@@ -6448 +6443 @@
                         if(values[i]>h_max){
                                 residual[i]=kappa*(values[i]-h_max);
-                                if(isefficientlayer) activate[i]=1.0;
                         }
                         else{
-                                if(isefficientlayer){
-                                        if(activeEpl[i]==1.0)activate[i]=1.0;
-                                }                               
                                 residual[i]=0.0;
                         }
                 }
         }
-        else{
-                for(int i=0;i<NUMVERTICES;i++){ 
-                        if(isefficientlayer){
-                                if(activeEpl[i]==1.0)activate[i]=1.0;
-                        }                               
-                }
-        }
 
         /*Add input to the element: */
-        if(isefficientlayer) this->inputs->AddInput(new TriaP1Input(HydrologydcMaskEplactiveEnum,&activate[0]));
         this->inputs->AddInput(new TriaP1Input(SedimentHeadEnum,values));
         this->inputs->AddInput(new TriaP1Input(SedimentHeadResidualEnum,residual));
@@ -6473 +6456 @@
         /*Free ressources:*/
         xDelete<int>(doflist);
-
 }
 /*}}}*/
@@ -6584 +6566 @@
 
         /*Constants*/
-        const int  numnodes         =NUMVERTICES;
-        int        flag;
+        const int  numnodes = NUMVERTICES;
+        IssmDouble flag     = 0.;
         IssmDouble active[numnodes];
         
         GetInputListOnVertices(&active[0],HydrologydcMaskEplactiveEnum);
-                
+
         for(int i=0;i<numnodes;i++) flag+=active[i];
         
         if(flag>0.){
-                for(int i=0;i<numnodes;i++) nodes[i]->Activate();
-        }
-        
-        else if(!this->AnyActive()){
-                for(int i=0;i<numnodes;i++) nodes[i]->Deactivate();
-                
+                for(int i=0;i<numnodes;i++){
+                        active_vec->SetValue(nodes[i]->Sid(),1.,INS_VAL);
+                }
         }
         else{
@@ -6617 +6596 @@
         IssmDouble  sedhead[numdof];
         IssmDouble  eplhead[numdof];
-
+        IssmDouble  residual[numdof];
 
         GetInputListOnVertices(&old_active[0],HydrologydcMaskEplactiveEnum);    
         GetInputListOnVertices(&sedhead[0],SedimentHeadEnum);
         GetInputListOnVertices(&eplhead[0],EplHeadEnum);
-
+        GetInputListOnVertices(&residual[0],SedimentHeadResidualEnum);
+
+        /*Get minimum sediment head*/
         sedheadmin=sedhead[0];
-        for(i=1;i<numdof;i++){
-                if(sedhead[i]<=sedheadmin)sedheadmin=sedhead[i];
-        }
+        for(i=1;i<numdof;i++) if(sedhead[i]<=sedheadmin)sedheadmin=sedhead[i];
         
         for(i=0;i<numdof;i++){
+
+                /*Activate EPL if residual is >0 */
+                if(residual[i]>0.){
+                        vec_mask->SetValue(nodes[i]->Sid(),1.,INS_VAL);
+                }
+                
                 /*If mask was alread one, keep one*/
-                if(old_active[i]>0.){
+                else if(old_active[i]>0.){
                         vec_mask->SetValue(nodes[i]->Sid(),1.,INS_VAL);
                 }
-                
-                this->GetHydrologyDCInefficientHmax(&h_max,nodes[i]);
+
                 /*Increase the size of the efficient system if needed*/
                 /*Increase is needed if the epl head reach the maximum value (sediment value for now)*/
+                this->GetHydrologyDCInefficientHmax(&h_max,nodes[i]);
                 if(eplhead[i]>=h_max){
-                        vec_mask->SetValue(nodes[i]->Sid(),1.,INS_VAL);
-                        for(j=0;j<numdof;j++){
-
-                                /*Increase of the domain is on the downstream node in term of sediment head*/
-                                if(sedhead[j] == sedheadmin){
-                                        vec_mask->SetValue(nodes[j]->Sid(),1.,INS_VAL);
-                                        break;
-                                }
-                        }
+                        //vec_mask->SetValue(nodes[i]->Sid(),1.,INS_VAL);
+                        //for(j=0;j<numdof;j++){
+
+                        //      /*Increase of the domain is on the downstream node in term of sediment head*/
+                        //      if(sedhead[j] == sedheadmin){
+                        //              vec_mask->SetValue(nodes[j]->Sid(),1.,INS_VAL);
+                        //              break;
+                        //      }
+                        //}
                 }
         }