KillIcebergsx.cpp

Go to the documentation of this file.

 
#include "./KillIcebergsx.h"
#include "../../shared/shared.h"
#include "../../toolkits/toolkits.h"
 
#include "../InputUpdateFromVectorx/InputUpdateFromVectorx.h"
 
void KillIcebergsx(FemModel* femmodel){
 
   /*Intermediaries*/
   int lid;
 
   /*retrieve vertex info and prepare element flag to speed up process*/
   int         nbv_local    = femmodel->vertices->Size();
   IssmDouble *local_mask   = xNewZeroInit<IssmDouble>(nbv_local);
   bool       *element_flag = xNewZeroInit<bool>(femmodel->elements->Size());
 
   /*Step 1, go through all elements and put 1 in local_mask if the element is grounded*/
   for(int i=0;i<femmodel->elements->Size();i++){
      Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
      if(!element->IsIceInElement()){
         /*Nothing to do, just flag element to speed up the computation*/
         element_flag[i] = true;
      }
      else{
         if(element->IsGrounded()){
            int numvertices = element->GetNumberOfVertices();
            for(int v=0;v<numvertices;v++) local_mask[element->vertices[v]->Lid()] = 1.;
         }
      }
   }
 
   /*Now we have 2 loops, one across cpus, and one for each cpus: we are going
    * to propagate the mask if an element is connected to a positive mask
    * already.  We then communicate to the other partitions. We stop when the
    * mask stops changing*/
   bool keepsyncing = true;
   while(keepsyncing){
 
      /*Get local mask from parallel vector*/
      femmodel->SyncLocalVectorWithClonesVerticesAdd(local_mask);
 
      /*Local iterations on partition*/
      bool keepgoing    = true;
      int  didsomething = 0;
      int  iter         = 1;
      while(keepgoing){
         //_printf0_("   -- Kill icebergs: local iteration "<<iter<<"\n");
 
         keepgoing    = false;
         didsomething = 0;
         for(int i=0;i<femmodel->elements->Size();i++){
            Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
 
            if(!element_flag[i]){
               int numvertices = element->GetNumberOfVertices();
               bool found1 = false;
               for(int j=0;j<numvertices;j++){
                  lid = element->vertices[j]->Lid();
                  if(local_mask[lid]>0.){
                     found1 = true;
                     break;
                  }
               }
               if(found1){
                  element_flag[i] = true;
                  for(int j=0;j<numvertices;j++){
                     lid = element->vertices[j]->Lid();
                     if(local_mask[lid]==0.){
                        local_mask[lid]=1.;
                        keepgoing = true;
                        didsomething = 1;
                     }
                  }
               }
            }
         }
         iter++;
      }
 
      /*Check how many iterations all cpus did*/
      int iter_max;
      ISSM_MPI_Reduce(&iter,&iter_max,1,ISSM_MPI_INT,ISSM_MPI_MAX,0,IssmComm::GetComm());
      ISSM_MPI_Bcast(&iter_max,1,ISSM_MPI_INT,0,IssmComm::GetComm());
      if(iter_max==2){
         /*If iter is only 2, nothing else was changed in the while loop above (iter is initialized as 1 and then ++)*/
         keepsyncing = false;
      }
   }
 
   /*Cleanup*/
   xDelete<bool>(element_flag);
 
   /*OK, now deactivate iceberg and count the number of deactivated vertices*/
   for(int i=0;i<femmodel->elements->Size();i++){
      Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
 
      if(element->IsIceInElement()){
         int  numvertices = element->GetNumberOfVertices();
         bool deactivate = false;
         for(int j=0;j<numvertices;j++){
            lid = element->vertices[j]->Lid();
            if(local_mask[lid]==0.){
               deactivate = true;
               break;
            }
         }
 
         if(deactivate){
            int  numvertices = element->GetNumberOfVertices();
            IssmDouble* values = xNew<IssmDouble>(numvertices);
            for(int j=0;j<numvertices;j++) values[j] = 1.; /*Anything >0 = no ice*/
            element->AddInput2(MaskIceLevelsetEnum,values,P1Enum);
            xDelete<IssmDouble>(values);
         }
      }
   }
 
   /*cleanup*/
   xDelete<IssmDouble>(local_mask);
}