ElementsAndVerticesPartitioning.cpp

Go to the documentation of this file.

 
#ifdef HAVE_CONFIG_H
   #include <config.h>
#else
#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
#endif
 
#include <string.h>
#include "../../classes/classes.h"
#include "../../shared/shared.h"
#include "../MeshPartitionx/MeshPartitionx.h"
#include "../ModelProcessorx/ModelProcessorx.h"
 
void  ElementsAndVerticesPartitioning(IoModel* iomodel){
 
   int numberofelements2d;
   int numberofvertices2d;
   int numlayers;
 
   /*intermediary: */
   int *epart          = NULL; //element partitioning.
   int *npart          = NULL; //node partitioning.
   int  elements_width;        //number of columns in elements (2d->3, 3d->6)
   int *elements2d     = NULL;
 
   /*Get my_rank:*/
   int my_rank   = IssmComm::GetRank();
   int num_procs = IssmComm::GetSize();
 
   /*First, check that partitioning has not yet been carryed out. Just check whether my_elements pointers is not already assigned a value: */
   if(iomodel->my_elements) return;
 
   /*Number of vertices per elements, needed to correctly retrieve data: */
   /*Determine parallel partitioning of elements: we use Metis for now. First load the data, then partition*/
   switch(iomodel->meshelementtype){
      case TriaEnum:
         elements_width=3;
         numberofelements2d = 0;
         numberofvertices2d = 0;
         numlayers          = 0;
         break;
      case TetraEnum:
         elements_width=4;
         numberofelements2d = 0;
         numberofvertices2d = 0;
         numlayers          = 0;
         break;
      case PentaEnum:
         elements_width=6;
         iomodel->FetchData(&elements2d,NULL,NULL,"md.mesh.elements2d");
         iomodel->FindConstant(&numberofelements2d,"md.mesh.numberofelements2d");
         iomodel->FindConstant(&numberofvertices2d,"md.mesh.numberofvertices2d");
         iomodel->FindConstant(&numlayers,"md.mesh.numberoflayers");
         break;
      default:
         _error_("mesh elements "<< EnumToStringx(iomodel->meshelementtype) <<" not supported yet");
   }
 
   /*Use ice levelset for weights*/
   int fordan = 0;
   int* weights = NULL;
   if(fordan){
      IssmDouble* icelevelset = NULL;
      iomodel->FetchData(&icelevelset,NULL,NULL,"md.mask.ice_levelset");
 
      weights = xNew<int>(iomodel->numberofvertices);
      for(int i=0;i<iomodel->numberofvertices;i++){
         if(icelevelset[i]>=0) weights[i] = 1;
         if(icelevelset[i]<0)  weights[i] = 100;
      }
      xDelete<IssmDouble>(icelevelset);
   }
 
   /*Partition and free resouces*/
   MeshPartitionx(&epart,&npart,iomodel->numberofelements,iomodel->numberofvertices,iomodel->elements,numberofelements2d,numberofvertices2d,elements2d,weights,numlayers,elements_width,iomodel->meshelementtype,num_procs);
   xDelete<int>(elements2d);
   xDelete<int>(npart);
   xDelete<int>(weights);
 
   if(fordan){
      for(int i=0;i<IssmComm::GetSize();i++){
         if(i==IssmComm::GetRank()){
            int temp =0;
            for(int j=0;j<iomodel->numberofelements;j++) if(epart[j]==i) temp++;
            _printf_("Partition #"<<i<<" number of elements: "<<temp<<"\n");
         }
         ISSM_MPI_Barrier(IssmComm::GetComm());
      }
   }
 
   /*Deal with rifts, they have to be included into one partition only, not several: */
   int numrifts;
   iomodel->FindConstant(&numrifts,"md.rifts.numrifts");
   if(numrifts){
      IssmDouble *riftinfo = NULL;
      iomodel->FetchData(&riftinfo,&numrifts,NULL,"md.rifts.riftstruct");
      for(int i=0;i<numrifts;i++){
         const int RIFTINFOSIZE = 12;
         int el1=reCast<int>(*(riftinfo+RIFTINFOSIZE*i+2))-1; //matlab indexing to c indexing
         int el2=reCast<int>(*(riftinfo+RIFTINFOSIZE*i+3))-1; //matlab indexing to c indexing
         epart[el2]=epart[el1]; //ensures that this pair of elements will be in the same partition, as well as the corresponding vertices;
      }
      iomodel->DeleteData(riftinfo,"md.rifts.riftstruct");
   }
 
   /*Create my_elements, used by each partition */
   bool *my_elements = xNewZeroInit<bool>(iomodel->numberofelements);
 
   /*Start figuring out, out of the partition, which elements belong to this cpu: */
   for(int i=0;i<iomodel->numberofelements;i++){
 
      if(my_rank==epart[i]) my_elements[i]=true;
   }
 
   /*Assign pointers to iomodel*/
   iomodel->epart      =epart;
   iomodel->my_elements=my_elements;
}