Charm++: Charm Source Code Documentation

00001 #include "ckvector3d.h"
00002 #include "charm-api.h"
00003 #include "refine.h"
00004 #include "fem.h"
00005 #include "fem_mesh.h"
00006 #include "femrefine.h"
00007 #include "ckhashtable.h"
00008 #include <assert.h>
00009 
00010 
00011 #define DEBUGINT(x) x
00012 //#define DEBUGINT(x) 
00013 
00014 
00015 class intdual{
00016 private:
00017   int x,y;
00018 public:
00019   intdual(int _x,int _y){
00020     if(_x <= _y){
00021       x = _x; y=_y;
00022     }else{
00023       x = _y; y= _x;
00024     }
00025   }
00026   inline int getx(){return x;};
00027   inline int gety(){return y;};
00028   inline CkHashCode hash() const {
00029     return (CkHashCode)(x+y);
00030   }
00031   static CkHashCode staticHash(const void *k,size_t){
00032     return ((intdual *)k)->hash();
00033   }
00034   inline int compare(intdual &t) const{
00035     return (t.getx() == x && t.gety() == y);
00036   }
00037   static int staticCompare(const void *a,const void *b,size_t){
00038     return ((intdual *)a)->compare((*(intdual *)b));
00039   }
00040 };
00041 
00042 void FEM_REFINE2D_Init(){
00043   REFINE2D_Init();  
00044 }
00045 
00046 FLINKAGE void FTN_NAME(FEM_REFINE2D_INIT,fem_refine2d_init)(void)
00047 {
00048   FEM_REFINE2D_Init();
00049 }
00050 
00051 
00052 void FEM_REFINE2D_Newmesh(int meshID,int nodeID,int elemID,int nodeBoundary){
00053   int nelems = FEM_Mesh_get_length(meshID,elemID);
00054   int nghost = FEM_Mesh_get_length(meshID,elemID+FEM_GHOST);
00055   int total = nghost + nelems;
00056   int *tempMesh = new int[3*total];
00057   int nnodes = FEM_Mesh_get_length(meshID,nodeID);
00058   int *tempBoundaries=NULL;
00059   if(nodeBoundary){
00060     tempBoundaries = new int[nnodes];
00061   }
00062   FEM_Mesh_data(meshID,elemID,FEM_CONN,&tempMesh[0],0,nelems,FEM_INDEX_0,3);
00063   FEM_Mesh_data(meshID,elemID+FEM_GHOST,FEM_CONN,&tempMesh[3*nelems],0,nghost,FEM_INDEX_0,3);
00064   
00065   for(int t=nelems;t<total;t++){
00066     for(int j=0;j<3;j++){
00067       if(FEM_Is_ghost_index(tempMesh[3*t+j])){
00068     tempMesh[3*t+j] += nelems;
00069       }
00070     }   
00071   }
00072   
00073   /*Set up the global ID's, for refinement*/
00074   int myID = FEM_My_partition();
00075   int *gid=new int[2*total];
00076   for (int i=0;i<nelems;i++) {
00077     gid[2*i+0]=myID; //Local element-- my chunk
00078     gid[2*i+1]=i; //Local number
00079   }
00080   int gid_fid=FEM_Create_field(FEM_INT,2,0,2*sizeof(int));
00081   FEM_Update_ghost_field(gid_fid,0,gid);
00082   
00083   if(nodeBoundary){
00084     FEM_Mesh_data(meshID,nodeID,FEM_BOUNDARY,tempBoundaries,0,nnodes,FEM_INT,1);
00085     printf("NODE BOUNDARIES-------------\n");
00086     for(int i=0;i<nnodes;i++){
00087       printf("%d %d \n",i,tempBoundaries[i]);
00088     }
00089     printf("------------------------\n");
00090   }
00091   
00092   /*Set up refinement framework*/
00093   /*FIX ME!  PASS IN EDGE BOUNDARIES! */
00094   const int *edgeConn = NULL;   
00095   const int *edgeBounds = NULL;
00096   int nEdges = 0;
00097   REFINE2D_NewMesh(meshID,nelems,total,nnodes,(int *)tempMesh,gid,tempBoundaries,edgeBounds, edgeConn, nEdges);
00098   if(tempBoundaries){
00099     delete [] tempBoundaries;
00100   }
00101   delete [] gid;
00102   delete [] tempMesh;
00103 }
00104 
00105 FLINKAGE void FTN_NAME(FEM_REFINE2D_NEWMESH,fem_refine2d_newmesh)(int *meshID,int *nodeID,int *elemID, int *nodeBoundary)
00106 {
00107   FEM_REFINE2D_Newmesh(*meshID,*nodeID,*elemID,*nodeBoundary);
00108 }
00109 
00110 class FEM_Refine_Operation_Data{
00111 public:
00112   int meshID,nodeID;
00113   int cur_nodes;
00114   CkVec<double> *coordVec;
00115   double *coord;
00116   FEM_Entity *node;
00117   CkVec<FEM_Attribute *> *attrs;
00118   FEM_Entity *elem;
00119   CkVec<FEM_Attribute *> *elemattrs;
00120   int sparseID,elemID;
00121   CkHashtableT<intdual,int> *newnodes;
00122   CkHashtableT<intdual,int> *nodes2sparse;
00123   FEM_Attribute *sparseConnAttr, *sparseBoundaryAttr;
00124   int nSplits;
00125   AllocTable2d<int> *connTable;
00126   FEM_Entity *sparse;
00127   CkVec<FEM_Attribute *> *sparseattrs;
00128   AllocTable2d<int> *validEdge;
00129   FEM_Refine_Operation_Data(){};
00130 };
00131 
00132 void FEM_Refine_Operation(FEM_Refine_Operation_Data *data,refineData &d);
00133 static int countValidEntities(int *validData,int total);
00134 
00135 void FEM_REFINE2D_Split(int meshID,int nodeID,double *coord,int elemID,double *desiredAreas,int sparseID){
00136   int nnodes = FEM_Mesh_get_length(meshID,nodeID);
00137   int nelems = FEM_Mesh_get_length(meshID,elemID);
00138   int actual_nodes = nnodes, actual_elems = nelems;
00139   
00140   FEM_Refine_Operation_Data refine_data;
00141   refine_data.meshID = meshID;
00142   refine_data.nodeID = nodeID;
00143   refine_data.sparseID = sparseID;
00144   refine_data.elemID = elemID;
00145   refine_data.cur_nodes = FEM_Mesh_get_length(meshID,nodeID);
00146   
00147   /*Copy the cordinates of the nodes into a vector, 
00148     the cordinates of the new nodes will be inserted
00149     into this vector and will be used to sort all the
00150     nodes on the basis of the distance from origin
00151   */
00152   CkVec<double> coordVec;
00153   for(int i=0;i<nnodes*2;i++){
00154     coordVec.push_back(coord[i]);
00155   }
00156   refine_data.coordVec = &coordVec;
00157   refine_data.coord = coord;
00158   
00159   /*find out the attributes of the node */
00160   FEM_Entity *e=refine_data.node = FEM_Entity_lookup(meshID,nodeID,"REFINE2D_Mesh");
00161   CkVec<FEM_Attribute *> *attrs = refine_data.attrs = e->getAttrVec();
00162   /* 
00163   FEM_DataAttribute *boundaryAttr = (FEM_DataAttribute *)e->lookup(FEM_BOUNDARY,"split");
00164   if(boundaryAttr != NULL){
00165     AllocTable2d<int> &boundaryTable = boundaryAttr->getInt();
00166     printf(" Node Boundary flags \n");
00167     for(int i=0;i<nnodes;i++){
00168         printf("Node %d flag %d \n",i,boundaryTable[i][0]);
00169     }
00170   }
00171   */
00172   FEM_Entity *elem = refine_data.elem = FEM_Entity_lookup(meshID,elemID,"REFIN2D_Mesh_elem");
00173   CkVec<FEM_Attribute *> *elemattrs = refine_data.elemattrs = elem->getAttrVec();
00174   
00175   FEM_Attribute *connAttr = elem->lookup(FEM_CONN,"REFINE2D_Mesh");
00176   if(connAttr == NULL){
00177     CkAbort("Grrrr element without connectivity \n");
00178   }
00179   AllocTable2d<int> &connTable = ((FEM_IndexAttribute *)connAttr)->get();
00180   refine_data.connTable = &connTable;
00181   
00182   //hashtable to store the new node number as a function of the two old numbers
00183   CkHashtableT<intdual,int> newnodes(nnodes);
00184   refine_data.newnodes = &newnodes;
00185   
00186   /* Get the FEM_BOUNDARY data of sparse elements and load it into a hashtable
00187     indexed by the 2 node ids that make up the edge. The data in the hashtable
00188     is the index number of the sparse element */
00189   FEM_Entity *sparse;
00190   CkVec<FEM_Attribute *> *sparseattrs;
00191   FEM_Attribute *sparseConnAttr, *sparseBoundaryAttr;
00192   AllocTable2d<int> *sparseConnTable, *sparseBoundaryTable;
00193   CkHashtableT<intdual,int> nodes2sparse(nelems);
00194   refine_data.nodes2sparse = &nodes2sparse;
00195   
00196   if(sparseID != -1){
00197     sparse = refine_data.sparse = FEM_Entity_lookup(meshID,sparseID,"REFINE2D_Mesh_sparse");
00198     refine_data.sparseattrs = sparseattrs = sparse->getAttrVec();
00199     refine_data.sparseConnAttr = sparseConnAttr = sparse->lookup(FEM_CONN,"REFINE2D_Mesh_sparse");
00200     sparseConnTable = &(((FEM_IndexAttribute *)sparseConnAttr)->get());
00201     refine_data.sparseBoundaryAttr = sparseBoundaryAttr = sparse->lookup(FEM_BOUNDARY,"REFINE2D_Mesh_sparse");
00202     if(sparseBoundaryAttr == NULL){
00203       CkAbort("Specified sparse elements without boundary conditions");
00204     }
00205     FEM_DataAttribute *validEdgeAttribute = (FEM_DataAttribute *)sparse->lookup(FEM_VALID,"REFINE2D_Mesh_sparse");
00206     if(validEdgeAttribute){
00207       refine_data.validEdge = &(validEdgeAttribute->getInt());
00208     }else{
00209       refine_data.validEdge = NULL;
00210     }
00211     /* since the default value in the hashtable is 0, to distinguish between 
00212        uninserted keys and the sparse element with index 0, the index of the 
00213        sparse elements is incremented by 1 while inserting. */
00214     //      printf("[%d] Sparse elements\n",FEM_My_partition());
00215     for(int j=0;j<sparse->size();j++){
00216       if(refine_data.validEdge == NULL || (*(refine_data.validEdge))[j][0]){
00217     int *cdata = (*sparseConnTable)[j];
00218     //      printf("%d < %d,%d > \n",j,cdata[0],cdata[1]);
00219     nodes2sparse.put(intdual(cdata[0],cdata[1])) = j+1;
00220       } 
00221     }
00222   }else{
00223     printf("Edge boundary conditions not passed into FEM_REFINE2D_Split \n");
00224   }
00225   
00226   //count the actual number of nodes and elements
00227   if(refine_data.node->lookup(FEM_VALID,"refine2D_splilt") != NULL){
00228     AllocTable2d<int> &validNodeTable = ((FEM_DataAttribute *)(refine_data.node->lookup(FEM_VALID,"refine2D_splilt")))->getInt();
00229     actual_nodes = countValidEntities(validNodeTable.getData(),nnodes);
00230   }
00231   if(refine_data.elem->lookup(FEM_VALID,"refine2D_splilt") != NULL){
00232     AllocTable2d<int> &validElemTable = ((FEM_DataAttribute *)(refine_data.elem->lookup(FEM_VALID,"refine2D_splilt")))->getInt();
00233     actual_elems = countValidEntities(validElemTable.getData(),nelems);
00234   }
00235   
00236   DEBUGINT(printf("%d %d \n",nnodes,nelems));   
00237   REFINE2D_Split(actual_nodes,coord,actual_elems,desiredAreas,&refine_data);
00238   
00239   int nSplits= refine_data.nSplits = REFINE2D_Get_Split_Length();
00240   DEBUGINT(printf("called REFINE2D_Split nSplits = %d \n",nSplits));
00241   
00242   if(nSplits == 0){
00243     return;
00244   }
00245   
00246   for(int split = 0;split < nSplits;split++){
00247     refineData op;
00248     REFINE2D_Get_Split(split,&op);
00249     FEM_Refine_Operation(&refine_data,op);
00250   }
00251   
00252   DEBUGINT(printf("Cordinate list length %d according to FEM %d\n",coordVec.size()/2,FEM_Mesh_get_length(meshID,nodeID)));
00253   IDXL_Sort_2d(FEM_Comm_shared(meshID,nodeID),coordVec.getVec());
00254   int read = FEM_Mesh_is_get(meshID) ;
00255   assert(read);
00256 }
00257 
00258 extern void splitEntity(IDXL_Side &c,
00259     int localIdx,int nBetween,int *between,int idxbase);
00260 
00261 void FEM_Modify_IDXL(FEM_Refine_Operation_Data *data,refineData &op){
00262   FEM_Node *node = (FEM_Node *)data->node;
00263   int between[2];
00264   between[0] = op.A;
00265   between[1] = op.B;
00266   splitEntity(node->shared,op.D,2,between,0);
00267 };
00268 
00269 void FEM_Refine_Operation(FEM_Refine_Operation_Data *data,refineData &op){
00270   
00271   int meshID = data->meshID;
00272   int nodeID = data->nodeID;
00273   int sparseID = data->sparseID;
00274   CkVec<FEM_Attribute *> *attrs = data->attrs;
00275   CkVec<FEM_Attribute *> *elemattrs = data->elemattrs;
00276   CkHashtableT<intdual,int> *newnodes=data->newnodes;
00277   CkHashtableT<intdual,int> *nodes2sparse = data->nodes2sparse;
00278   FEM_Attribute *sparseConnAttr = data->sparseConnAttr, *sparseBoundaryAttr = data->sparseBoundaryAttr;
00279   AllocTable2d<int> *connTable = data->connTable;
00280   double *coord = data->coord;
00281   AllocTable2d<int> *sparseConnTable, *sparseBoundaryTable;
00282   CkVec<FEM_Attribute *> *sparseattrs = data->sparseattrs;
00283   /*
00284   FEM_DataAttribute *boundaryAttr = (FEM_DataAttribute *)data->node->lookup(FEM_BOUNDARY,"split");
00285   if(boundaryAttr != NULL){
00286   AllocTable2d<int> &boundaryTable = boundaryAttr->getInt();
00287   printf(" Node Boundary flags \n");
00288   for(int i=0;i<data->node->size();i++){
00289   printf("Node %d flag %d \n",i,boundaryTable[i][0]);
00290   }
00291   }
00292   */
00293   
00294   int tri=op.tri,A=op.A,B=op.B,C=op.C,D=op.D;
00295   double frac=op.frac;
00296   // current number of nodes in the mesh
00297   int *connData = connTable->getData();
00298   int flags=op.flag;
00299     
00300   if((flags & 0x1) || (flags & 0x2)){
00301     //new node 
00302     DEBUGINT(CkPrintf("---- Adding node %d\n",D));          
00303     /*  lastA=A; lastB=B; lastD=D;*/
00304     if (A>=data->cur_nodes) CkAbort("Calculated A is invalid!");
00305     if (B>=data->cur_nodes) CkAbort("Calculated B is invalid!");
00306     if(D >= data->cur_nodes){
00307       data->node->setLength(D+1);
00308       data->cur_nodes = D+1;
00309     }   
00310     for(int i=0;i<attrs->size();i++){
00311       FEM_Attribute *a = (FEM_Attribute *)(*attrs)[i];
00312       if(a->getAttr()<FEM_ATTRIB_TAG_MAX){
00313     FEM_DataAttribute *d = (FEM_DataAttribute *)a;
00314     d->interpolate(A,B,D,frac);
00315       }else{
00316     /* Boundary value of new node D should be boundary value of the edge
00317        (sparse element) that contains the two nodes A & B. */
00318     if(a->getAttr() == FEM_BOUNDARY){
00319       if(sparseID != -1){
00320         int sidx = nodes2sparse->get(intdual(A,B))-1;
00321         if(sidx == -1){
00322           CkAbort("no sparse element between these 2 nodes, are they really connected ??");
00323         }
00324         sparseBoundaryTable = &(((FEM_DataAttribute *)sparseBoundaryAttr)->getInt());
00325         
00326         int boundaryVal = ((*sparseBoundaryTable)[sidx])[0];
00327         
00328         (((FEM_DataAttribute *)a)->getInt()[D])[0] = boundaryVal;
00329       }else{
00330         /* if sparse elements don't exist, just do simple interpolation */
00331         FEM_DataAttribute *d = (FEM_DataAttribute *)a;
00332         d->interpolate(A,B,D,frac);
00333       }
00334     }
00335       } 
00336     }
00337     int AandB[2];
00338     AandB[0]=A;
00339     AandB[1]=B;
00340     /* Add a new node D between A and B */
00341     IDXL_Add_entity(FEM_Comm_shared(meshID,nodeID),D,2,AandB);
00342     double Dx = coord[2*A]*(1-frac)+frac*coord[2*B];
00343     double Dy = coord[2*A+1]*(1-frac)+frac*coord[2*B+1];
00344     data->coordVec->push_back(Dx);
00345     data->coordVec->push_back(Dy);
00346     newnodes->put(intdual(A,B))=D;
00347     /* add the new sparse element <D,B> and modify the connectivity of the old
00348        one from <A,B> to <A,D> and change the hashtable to reflect that change
00349     */
00350     if(sparseID != -1){
00351       int oldsidx = nodes2sparse->get(intdual(A,B))-1;
00352       int newsidx = data->sparse->size();
00353       data->sparse->setLength(newsidx+1);
00354       for(int satt = 0;satt<sparseattrs->size();satt++){
00355     if((*sparseattrs)[satt]->getAttr() == FEM_CONN){
00356       /* change the conn of the old sparse to A,D and new one to B,D */
00357       sparseConnTable = &(((FEM_IndexAttribute *)sparseConnAttr)->get());
00358       int *oldconn = (*sparseConnTable)[oldsidx];
00359       int *newconn = (*sparseConnTable)[newsidx];
00360       oldconn[0] = A;
00361       oldconn[1] = D;
00362       newconn[0] = D;
00363       newconn[1] = B;
00364       //printf("<%d,%d> edge being split into <%d,%d> <%d,%d> \n",A,B,A,D,D,B);
00365     }else{
00366       /* apart from conn copy everything else */
00367       FEM_Attribute *attr = (FEM_Attribute *)(*sparseattrs)[satt];
00368       attr->copyEntity(newsidx,*attr,oldsidx);
00369     }
00370       }
00371       /* modify the hashtable - delete the old edge and the new ones */
00372       nodes2sparse->remove(intdual(A,B));
00373       nodes2sparse->put(intdual(A,D)) = oldsidx+1;
00374       nodes2sparse->put(intdual(D,B)) = newsidx+1;
00375     }
00376   }
00377   //add a new triangle
00378   /*TODO: replace  FEM_ELEM with parameter*/
00379   int newTri =  op._new;
00380   int cur_elements  = FEM_Mesh_get_length(data->meshID,data->elemID);
00381   DEBUGINT(CkPrintf("---- Adding triangle %d after splitting %d \n",newTri,tri));
00382   if(newTri >= cur_elements){
00383     data->elem->setLength(newTri+1);
00384   } 
00385   D = newnodes->get(intdual(A,B));
00386   
00387   for(int j=0;j<elemattrs->size();j++){
00388     if((*elemattrs)[j]->getAttr() == FEM_CONN){
00389       DEBUGINT(CkPrintf("elem attr conn code %d \n",(*elemattrs)[j]->getAttr()));
00390       //it is a connectivity attribute.. get the connectivity right
00391       FEM_IndexAttribute *connAttr = (FEM_IndexAttribute *)(*elemattrs)[j];
00392       AllocTable2d<int> &table = connAttr->get();
00393       DEBUGINT(CkPrintf("Table of connectivity attribute starts at %p width %d \n",table[0],connAttr->getWidth()));
00394       int *oldRow = table[tri];
00395       int *newRow = table[newTri];
00396       for (int i=0;i<3;i++){
00397     if (oldRow[i]==A){
00398       oldRow[i]=D;  
00399       DEBUGINT(CkPrintf("In triangle %d %d replaced by %d \n",tri,A,D));
00400     }
00401       } 
00402       for (int i=0; i<3; i++) {
00403     if (oldRow[i] == B){
00404       newRow[i] = D;
00405     }   
00406     else if (oldRow[i] == C){
00407       newRow[i] = C;
00408     }   
00409     else if (oldRow[i] == D){
00410       newRow[i] = A;
00411     }   
00412       }
00413       DEBUGINT(CkPrintf("New Triangle %d  (%d %d %d) conn %p\n",newTri,newRow[0],newRow[1],newRow[2],newRow));
00414     }else{
00415       FEM_Attribute *elattr = (FEM_Attribute *)(*elemattrs)[j];
00416       if(elattr->getAttr() < FEM_ATTRIB_FIRST){ 
00417     elattr->copyEntity(newTri,*elattr,tri);
00418       }
00419     }
00420   }
00421   
00422   if(sparseID != -1){ /* add sparse element (edge between C and D) */
00423     int cdidx = data->sparse->size();
00424     data->sparse->setLength(cdidx+1);
00425     for(int satt = 0; satt < sparseattrs->size();satt++){
00426       if((*sparseattrs)[satt]->getAttr() == FEM_CONN){
00427     sparseConnTable = &(((FEM_IndexAttribute *)sparseConnAttr)->get());
00428     int *cdconn = (*sparseConnTable)[cdidx];
00429     cdconn[0]=C;
00430     cdconn[1]=D;
00431       }
00432       if((*sparseattrs)[satt]->getAttr() == FEM_BOUNDARY){
00433     /* An edge connecting C and D has to be an internal edge */
00434     sparseBoundaryTable = &(((FEM_DataAttribute *)sparseBoundaryAttr)->getInt());
00435     ((*sparseBoundaryTable)[cdidx])[0] = 0;
00436       }
00437     }
00438     if(data->validEdge){
00439       (*(data->validEdge))[cdidx][0] = 1;
00440     }
00441     nodes2sparse->put(intdual(C,D)) = cdidx+1;
00442   }
00443 }
00444 
00445 
00446 FLINKAGE void FTN_NAME(FEM_REFINE2D_SPLIT,fem_refine2d_split)(int *meshID,int *nodeID,double *coord,int *elemID,double *desiredAreas){
00447     FEM_REFINE2D_Split(*meshID,*nodeID,coord,*elemID,desiredAreas);
00448 }
00449 
00450 FLINKAGE void FTN_NAME(FEM_REFINE2D_SPLIT_EDGE,fem_refine2d_split_edge)(int *meshID,int *nodeID,double *coord,int *elemID,double *desiredAreas,int *sparseID){
00451     FEM_REFINE2D_Split(*meshID,*nodeID,coord,*elemID,desiredAreas,*sparseID);
00452 }
00453 
00454 FLINKAGE void FTN_NAME(CMIMEMORYCHECK,cmimemorycheck)(){
00455     CmiMemoryCheck();
00456 }
00457 
00458 class FEM_Operation_Data{
00459 public:
00460   double *coord;
00461   int *connData,*nodeBoundaryData;
00462   int *validNodeData,*validElemData;
00463   FEM_Node *node;
00464   CkHashtableT<intdual,int> *nodes2sparse;
00465   AllocTable2d<int> *validEdge;
00466   AllocTable2d<int> *sparseConnTable;
00467   FEM_Attribute *sparseBoundaryAttr;
00468   FEM_Operation_Data(){
00469   };
00470 };
00471 void FEM_Coarsen_Operation(FEM_Operation_Data *coarsen_data, coarsenData &operation);
00472 
00473 void FEM_REFINE2D_Coarsen(int meshID,int nodeID,double *coord,int elemID,double *desiredAreas,int sparseID){
00474   int nnodes = FEM_Mesh_get_length(meshID,nodeID);
00475   int nelems = FEM_Mesh_get_length(meshID,elemID);
00476   int nodeCount=0,elemCount=0;
00477   
00478   /* The attributes of the different entities */
00479   FEM_Entity *node=FEM_Entity_lookup(meshID,nodeID,"REFINE2D_Mesh");
00480   CkVec<FEM_Attribute *> *attrs = node->getAttrVec();
00481   FEM_Attribute *validNodeAttr = node->lookup(FEM_VALID,"FEM_COARSEN");
00482   FEM_Attribute *nodeBoundaryAttr = node->lookup(FEM_BOUNDARY,"FEM_COARSEN");
00483   if(!nodeBoundaryAttr){
00484     printf("Warning:- no boundary flags for nodes \n");
00485   }
00486   
00487   FEM_Entity *elem = FEM_Entity_lookup(meshID,elemID,"REFIN2D_Mesh_elem");
00488   CkVec<FEM_Attribute *> *elemattrs = elem->getAttrVec();
00489   FEM_Attribute *validElemAttr = elem->lookup(FEM_VALID,"FEM_COARSEN");
00490   
00491   FEM_Attribute *connAttr = elem->lookup(FEM_CONN,"REFINE2D_Mesh");
00492   if(connAttr == NULL){
00493     CkAbort("Grrrr element without connectivity \n");
00494   }
00495   AllocTable2d<int> &connTable = ((FEM_IndexAttribute *)connAttr)->get();
00496   AllocTable2d<int>&validNodeTable = ((FEM_DataAttribute *)validNodeAttr)->getInt();
00497   AllocTable2d<int>&validElemTable = ((FEM_DataAttribute *)validElemAttr)->getInt();
00498   FEM_Operation_Data *coarsen_data = new FEM_Operation_Data;
00499   coarsen_data->node = (FEM_Node *)node;
00500   coarsen_data->coord = coord;
00501   int *connData = coarsen_data->connData= connTable.getData();
00502   int *validNodeData = coarsen_data->validNodeData = validNodeTable.getData();
00503   int *validElemData = coarsen_data->validElemData = validElemTable.getData();
00504   /* Extract the data for node boundaries */
00505   int *nodeBoundaryData= coarsen_data->nodeBoundaryData =NULL;
00506   if(nodeBoundaryAttr){
00507     AllocTable2d<int> &nodeBoundaryTable = ((FEM_DataAttribute *)nodeBoundaryAttr)->getInt();
00508     nodeBoundaryData = coarsen_data->nodeBoundaryData = nodeBoundaryTable.getData();
00509   }
00510   
00511   for(int k=0;k<nnodes;k++){
00512     int valid = validNodeData[k];
00513     if(validNodeData[k]){
00514       nodeCount++;
00515     }
00516   }
00517   for(int k=0;k<nelems;k++){
00518     if(validElemData[k]){
00519       elemCount++;
00520     }
00521   }
00522   
00523   /* populate the hashtable from nodes2edges with the valid edges */
00524   CkHashtableT<intdual,int> nodes2sparse(nelems);
00525   coarsen_data->nodes2sparse = &nodes2sparse;
00526   FEM_Attribute *sparseBoundaryAttr;
00527   if(sparseID != -1){
00528     FEM_Entity *sparse = FEM_Entity_lookup(meshID,sparseID,"Coarsen_sparse");
00529     FEM_DataAttribute *validEdgeAttribute = (FEM_DataAttribute *)sparse->lookup(FEM_VALID,"Coarsen_sparse");
00530     FEM_IndexAttribute *sparseConnAttribute = (FEM_IndexAttribute *)sparse->lookup(FEM_CONN,"Coarsen_sparse");
00531     AllocTable2d<int> *sparseConnTable = coarsen_data->sparseConnTable = &(sparseConnAttribute->get());
00532     coarsen_data->sparseBoundaryAttr = sparseBoundaryAttr = sparse->lookup(FEM_BOUNDARY,"REFINE2D_Mesh_sparse");
00533     if(sparseBoundaryAttr == NULL){
00534       CkAbort("Specified sparse elements without boundary conditions");
00535     }
00536 
00537     
00538     if(validEdgeAttribute){
00539       coarsen_data->validEdge = &(validEdgeAttribute->getInt());
00540     }else{
00541       coarsen_data->validEdge = NULL;
00542     }
00543     /* since the default value in the hashtable is 0, to distinguish between 
00544        uninserted keys and the sparse element with index 0, the index of the 
00545        sparse elements is incremented by 1 while inserting. */
00546     printf("[%d] Sparse elements\n",FEM_My_partition());
00547     for(int j=0;j<sparse->size();j++){
00548       if(coarsen_data->validEdge == NULL || (*(coarsen_data->validEdge))[j][0]){
00549     int *cdata = (*sparseConnTable)[j];
00550     printf("%d < %d,%d > \n",j,cdata[0],cdata[1]);
00551     nodes2sparse.put(intdual(cdata[0],cdata[1])) = j+1;
00552       } 
00553     }
00554   }else{
00555     coarsen_data->validEdge = NULL;
00556   }
00557   
00558   DEBUGINT(printf("coarsen %d %d \n",nodeCount,elemCount)); 
00559   REFINE2D_Coarsen(nodeCount,coord,elemCount,desiredAreas,coarsen_data);
00560   int nCollapses = REFINE2D_Get_Collapse_Length();
00561   
00562   
00563   for(int i=0;i<nnodes;i++){
00564     int valid = validNodeData[i];
00565   }
00566   delete coarsen_data;
00567 }  
00568 
00569 void interpolateNode(FEM_Node *node,int A,int B,int D,double frac);
00570 
00571 void FEM_Coarsen_Operation(FEM_Operation_Data *coarsen_data, coarsenData &operation){
00572   double *coord = coarsen_data->coord;
00573   int tri,nodeToThrow,nodeToKeep,n1,n2;
00574   int *connData = coarsen_data->connData;
00575   int *validNodeData = coarsen_data->validNodeData;
00576   int *validElemData = coarsen_data->validElemData;
00577   int *nodeBoundaryData = coarsen_data->nodeBoundaryData;
00578   FEM_Attribute *sparseBoundaryAttr = coarsen_data->sparseBoundaryAttr;
00579   AllocTable2d<int> *sparseBoundaryTable;
00580 
00581   switch(operation.type){
00582   case COLLAPSE: 
00583     {
00584       tri = operation.data.cdata.elemID;
00585       nodeToKeep = operation.data.cdata.nodeToKeep;
00586       nodeToThrow = operation.data.cdata.nodeToDelete;
00587       int opNode = 0;
00588       for (int i=0; i<3; i++) {
00589     if ((connData[3*tri+i] != nodeToThrow) &&
00590         (connData[3*tri+i] != nodeToKeep)) {
00591       opNode = connData[3*tri+i];
00592       break;
00593     }
00594       }
00595       CkPrintf("Collapse %d, nodeToKeep %d, nodeToThrow %d, opNode %d\n",
00596            tri, nodeToKeep, nodeToThrow, opNode);
00597       sparseBoundaryTable = &(((FEM_DataAttribute *)sparseBoundaryAttr)->getInt());
00598       int delEdgeIdx = coarsen_data->nodes2sparse->get(intdual(nodeToThrow,opNode))-1;
00599       int keepEdgeIdx = coarsen_data->nodes2sparse->get(intdual(nodeToKeep,opNode))-1;
00600       int delBC = ((*sparseBoundaryTable)[delEdgeIdx])[0];
00601       int keepBC = ((*sparseBoundaryTable)[keepEdgeIdx])[0];
00602       if (delBC > keepBC) {
00603     ((*sparseBoundaryTable)[keepEdgeIdx])[0] = delBC;
00604       }
00605       
00606       if(operation.data.cdata.flag & 0x1 || operation.data.cdata.flag & 0x2){
00607     interpolateNode(coarsen_data->node,nodeToKeep,nodeToThrow,nodeToKeep,operation.data.cdata.frac);
00608     coord[2*nodeToKeep] = operation.data.cdata.newX;
00609     coord[2*nodeToKeep+1] = operation.data.cdata.newY;
00610     validNodeData[nodeToThrow] = 0;
00611     validNodeData[nodeToKeep] = 1;
00612     DEBUGINT(printf("---------Collapse <%d,%d> invalidating node %d and element %d \n",nodeToKeep,nodeToThrow,nodeToThrow,tri));
00613     if(coarsen_data->validEdge){    
00614       int sidx = coarsen_data->nodes2sparse->get(intdual(nodeToKeep,nodeToThrow))-1;
00615       coarsen_data->nodes2sparse->remove(intdual(nodeToKeep,nodeToThrow));
00616       (*(coarsen_data->validEdge))[sidx][0] = 0;
00617       DEBUGINT(printf("---- Deleting edge %d between nodes %d and %d \n",sidx,nodeToKeep,nodeToThrow));
00618       if (delEdgeIdx >=0) {
00619         coarsen_data->nodes2sparse->remove(intdual(opNode,nodeToThrow));
00620         (*(coarsen_data->validEdge))[delEdgeIdx][0] = 0;
00621         DEBUGINT(printf("---- Deleting edge %d between nodes %d and %d \n",delEdgeIdx,opNode,nodeToThrow));
00622       }
00623     }   
00624       }
00625       else {
00626     if (delEdgeIdx >=0) {
00627       coarsen_data->nodes2sparse->remove(intdual(opNode,nodeToThrow));
00628       (*(coarsen_data->validEdge))[delEdgeIdx][0] = 0;
00629       DEBUGINT(printf("---- Deleting edge %d between nodes %d and %d \n",delEdgeIdx,opNode,nodeToThrow));
00630     }   
00631       }
00632       validElemData[tri] = 0;
00633       connData[3*tri] = connData[3*tri+1] = connData[3*tri+2] = -1;
00634     }
00635     break;
00636   case  UPDATE:
00637     if(validNodeData[operation.data.udata.nodeID]){
00638       coord[2*(operation.data.udata.nodeID)] = operation.data.udata.newX;
00639       coord[2*(operation.data.udata.nodeID)+1] = operation.data.udata.newY;
00640       if(nodeBoundaryData){
00641     nodeBoundaryData[operation.data.udata.nodeID]=operation.data.udata.boundaryFlag;
00642       }
00643     }else{
00644       DEBUGINT(printf("[%d] WEIRD -- update operation for invalid node %d \n",CkMyPe(),operation.data.udata.nodeID));
00645     }
00646     break;
00647   case REPLACE:
00648     if(validElemData[operation.data.rddata.elemID]){
00649       if(connData[3*operation.data.rddata.elemID+operation.data.rddata.relnodeID] == operation.data.rddata.oldNodeID){
00650     connData[3*operation.data.rddata.elemID+operation.data.rddata.relnodeID] = operation.data.rddata.newNodeID;
00651     if(validNodeData[operation.data.rddata.oldNodeID]){
00652       validNodeData[operation.data.rddata.oldNodeID]=0;
00653     }
00654     if(coarsen_data->validEdge){
00655       //remove the edges containing oldNodeID and add edges containing newNodeID in the nodes2sparse hashtable
00656       //update the connectivity information of the edges
00657       for(int i=0;i<3;i++){
00658         //find the two nodes apart from the one being replaced
00659         if(i != operation.data.rddata.relnodeID){
00660           int otherNode = connData[3*operation.data.rddata.elemID+i];
00661           int edgeIdx = coarsen_data->nodes2sparse->get(intdual(operation.data.rddata.oldNodeID,otherNode))-1;
00662           if(edgeIdx >= 0){
00663         //The edge connectivity has not been updated on this processor
00664         coarsen_data->nodes2sparse->remove(intdual(operation.data.rddata.oldNodeID,otherNode));
00665         DEBUGINT(printf("---- Deleting edge %d between nodes %d and %d \n",edgeIdx,operation.data.rddata.oldNodeID,otherNode));
00666         coarsen_data->nodes2sparse->put(intdual(operation.data.rddata.newNodeID,otherNode)) = edgeIdx+1;
00667         DEBUGINT(printf("---- Adding edge %d between nodes %d and %d \n",edgeIdx,operation.data.rddata.newNodeID,otherNode));
00668         (*(coarsen_data->sparseConnTable))[edgeIdx][0] = otherNode;
00669         (*(coarsen_data->sparseConnTable))[edgeIdx][1] = operation.data.rddata.newNodeID;
00670           } 
00671         }
00672       }
00673     }
00674       }else{
00675     DEBUGINT(printf("[%d] WEIRD -- REPLACE operation for element %d specifies different node number %d \n",CkMyPe(),operation.data.rddata.elemID,operation.data.rddata.oldNodeID));
00676       }
00677     }else{
00678       DEBUGINT(printf("[%d] WEIRD -- REPLACE operation for invalid element %d \n",CkMyPe(),operation.data.rddata.elemID));
00679     }
00680     DEBUGINT(printf("---------Replace invalidating node %d with %d in element %d\n",operation.data.rddata.oldNodeID,operation.data.rddata.newNodeID,operation.data.rddata.elemID));
00681     break;
00682   default:
00683     DEBUGINT(printf("[%d] WEIRD -- COARSENDATA type == invalid \n",CkMyPe()));
00684     CmiAbort("COARSENDATA type == invalid");
00685   }         
00686 };
00687 
00688 
00689 FLINKAGE void FTN_NAME(FEM_REFINE2D_COARSEN,fem_refine2d_coarsen)(int *meshID,int *nodeID,double *coord,int *elemID,double *desiredAreas,int *sparseID){
00690   FEM_REFINE2D_Coarsen(*meshID,*nodeID,coord,*elemID,desiredAreas,*sparseID);
00691 }
00692 
00693 void interpolateNode(FEM_Node *node,int A,int B,int D,double frac){
00694   CkVec<FEM_Attribute *>*attrs = node->getAttrVec();
00695   for(int i=0;i<attrs->size();i++){
00696     FEM_Attribute *a = (FEM_Attribute *)(*attrs)[i];
00697     if(a->getAttr()<FEM_ATTRIB_TAG_MAX){
00698       FEM_DataAttribute *d = (FEM_DataAttribute *)a;
00699       AllocTable2d<double> *doubleData;
00700       if(a->getAttr() == FEM_DATA+6){
00701     doubleData = &d->getDouble();
00702     printf("A %d %.6lf %.6lf B %d %.6lf %.6lf frac %.6lf ",A,(*doubleData)[A][0],(*doubleData)[A][1],B,(*doubleData)[B][0],(*doubleData)[B][1],frac);
00703       }
00704       d->interpolate(A,B,D,frac);
00705       if(a->getAttr() == FEM_DATA+6){
00706     printf("D %d %.6lf %.6lf \n",D,(*doubleData)[D][0],(*doubleData)[D][1]);
00707       }
00708     }   
00709   } 
00710 }
00711 
00712 
00713 static int countValidEntities(int *validData,int total){
00714   int sum =0 ;
00715   for(int i=0;i<total;i++){
00716     sum += validData[i];
00717   }
00718   return sum;
00719 }
libs/ck-libs/TMRC2D/femrefine.C
