Changeset 16788

Timestamp:

11/15/13 15:04:16 (11 years ago)

Author:

bdef

Message:

BUG+CHG: fixing a conditional jmp in pengrid and changing the way in which we deal with the active nodes for L2ProjectionEPL

Location:

issm/trunk-jpl/src/c/classes

Files:

• Elements/Penta.cpp (modified) (1 diff)
• Elements/Tria.cpp (modified) (2 diffs)
• FemModel.cpp (modified) (1 diff)
• Loads/Pengrid.cpp (modified) (2 diffs)

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

@@ -9519 +9519 @@
 void Penta::HydrologyEPLGetActive(Vector<IssmDouble>* active_vec){
 
-        if (!IsOnBed())return;
-
+        if (!IsOnBed()){
+                return;
+        }
         Tria* tria=(Tria*)SpawnTria(0); //lower face is 0, upper face is 1.
         tria->HydrologyEPLGetActive(active_vec);

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

@@ -6989 +6989 @@
                 
                 for(int i=0;i<numdof;i++){
-                        /*Keeping thickness to 1 if EPL is not active*/
-                        if(activeEpl[i]==0.0){
-                                thickness[i]=init_thick;
-                        }
-                        else{
-
+                        /*Keeping thickness to initial value if EPL is not active*/
+                        /* if(activeEpl[i]==0.0){ */
+                        /*      thickness[i]=init_thick; */
+                        /* } */
+                        /* else{ */
                                 /*Compute first the effective pressure in the EPL*/
                                 EPL_N=gravity*((rho_ice*ice_thickness[i])-(rho_water*(eplhead[i]-bed[i])));
@@ -7003 +7002 @@
                                 /*And proceed to the real thing*/
                                 thickness[i] = old_thickness[i]*(1+((rho_water*gravity*dt)/(rho_ice*latentheat))*epl_conductivity*pow(EPLgrad,2.0)-2.0*(A*dt/(pow(n,n)))*(pow(EPL_N,n)));
-                        }
+                                //}
                 }
                 this->inputs->AddInput(new TriaInput(HydrologydcEplThicknessEnum,thickness,P1Enum));

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

@@ -1423 +1423 @@
 void FemModel::UpdateConstraintsL2ProjectionEPLx(void){ /*{{{*/
 
-        for(int i=0;i<elements->Size();i++){
+        Vector<IssmDouble>* active        = NULL;
+        IssmDouble*         serial_active = NULL;
+
+        /*update node activity. If one element is connected to mask=1, all nodes are active*/
+        active=new Vector<IssmDouble>(nodes->NumberOfNodes(HydrologyDCEfficientAnalysisEnum));
+        for (int i=0;i<elements->Size();i++){
                 Element* element=dynamic_cast<Element*>(elements->GetObjectByOffset(i));
-                element->UpdateConstraintsL2ProjectionEPL();
-        }
+                element->HydrologyEPLGetActive(active);
+        }
+
+        /*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=dynamic_cast<Node*>(nodes->GetObjectByOffset(i));
+                if(node->InAnalysis(L2ProjectionEPLAnalysisEnum)){
+                        if(serial_active[node->Sid()]==1.){
+                                node->Activate();
+                                counter++;
+                        }
+                        else{
+                                node->Deactivate();
+                        }
+                }
+        }
+        xDelete<IssmDouble>(serial_active);
+        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;
+        if(VerboseSolution()) _printf0_("   Number of active nodes L2 Projection: "<< counter <<"\n");
+
+        /*Update dof indexings*/
+        this->UpdateConstraintsx();
+
+        /* for(int i=0;i<elements->Size();i++){ */
+        /*      Element* element=dynamic_cast<Element*>(elements->GetObjectByOffset(i)); */
+        /*      element->UpdateConstraintsL2ProjectionEPL(); */
+        /* } */
 
 }

issm/trunk-jpl/src/c/classes/Loads/Pengrid.cpp

@@ -648 +648 @@
                 return;
         }
+        if(!element->IsOnBed()){
+                unstable=0;
+                active=0;
+                *punstable=unstable;
+                return;
+        }
 
         /*Get sediment water head h*/
@@ -654 +660 @@
         parameters->FindParam(&penalty_lock,HydrologydcPenaltyLockEnum);
 
-        if (h>h_max)
+        if (h>h_max){
          new_active=1;
-        else
+        }
+        else{
          new_active=0;
+        }
 
         if(this->active==new_active){