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 
00008 class intdual{
00009         private:
00010                 int x,y;
00011         public:
00012                 intdual(int _x,int _y){
00013                         if(_x <= _y){
00014                                 x = _x; y=_y;
00015                         }else{
00016                                 x = _y; y= _x;
00017                         }
00018                 }
00019                 int getx(){return x;};
00020                 int gety(){return y;};
00021                 inline CkHashCode hash() const {
00022                         return (CkHashCode)(x+y);
00023                 }
00024                 static CkHashCode staticHash(const void *k,size_t){
00025                         return ((intdual *)k)->hash();
00026                 }
00027                 inline int compare(intdual &t) const{
00028                         return (t.getx() == x && t.gety() == y);
00029                 }
00030                 static int staticCompare(const void *a,const void *b,size_t){
00031                         return ((intdual *)a)->compare((*(intdual *)b));
00032                 }
00033 };
00034 
00035 void FEM_REFINE2D_Init(){
00036   REFINE2D_Init();      
00037 }
00038 
00039 FDECL void FTN_NAME(FEM_REFINE2D_INIT,fem_refine2d_init)(void)
00040 {
00041   FEM_REFINE2D_Init();
00042 }
00043 
00044 
00045 
00046 
00047 void FEM_REFINE2D_Newmesh(int meshID,int nodeID,int elemID){
00048         int nelems = FEM_Mesh_get_length(meshID,elemID);
00049         int     nghost = FEM_Mesh_get_length(meshID,elemID+FEM_GHOST);
00050         int total = nghost + nelems;
00051         int *tempMesh = new int[3*total];
00052         FEM_Mesh_data(meshID,elemID,FEM_CONN,&tempMesh[0],0,nelems,FEM_INDEX_0,3);
00053         FEM_Mesh_data(meshID,elemID+FEM_GHOST,FEM_CONN,&tempMesh[3*nelems],0,nghost,FEM_INDEX_0,3);
00054 
00055         for(int t=nelems;t<total;t++){
00056                 for(int j=0;j<3;j++){
00057                         if(FEM_Is_ghost_index(tempMesh[3*t+j])){
00058                                 tempMesh[3*t+j] += nelems;
00059                         }
00060                 }       
00061         }
00062         int maxnode=0,maxid=0;
00063         for(int i=0;i<3*total;i++){
00064                 if(tempMesh[i] > maxnode){
00065                         maxnode = tempMesh[i];
00066                         maxid=i;
00067                 }
00068         }
00069   /*Set up the global ID's, for refinement*/
00070         int myID = FEM_My_partition();
00071   int *gid=new int[2*total];
00072   for (int i=0;i<nelems;i++) {
00073     gid[2*i+0]=myID; //Local element-- my chunk
00074     gid[2*i+1]=i; //Local number
00075   }
00076   int gid_fid=FEM_Create_field(FEM_INT,2,0,2*sizeof(int));
00077   FEM_Update_ghost_field(gid_fid,0,gid);
00078         
00079   /*Set up refinement framework*/
00080   printf("NewMesh %d %d %d maxid %d \n",nelems,total,maxnode,maxid);
00081   REFINE2D_NewMesh(nelems,total,(int *)tempMesh,gid);
00082         delete [] gid;
00083         delete [] tempMesh;
00084 }
00085 
00086 FDECL void FTN_NAME(FEM_REFINE2D_NEWMESH,fem_refine2d_newmesh)(int *meshID,int *nodeID,int *elemID){
00087         FEM_REFINE2D_Newmesh(*meshID,*nodeID,*elemID);
00088 }
00089 
00090 
00091 
00092 void FEM_REFINE2D_Split(int meshID,int nodeID,double *coord,int elemID,double *desiredAreas,int sparseID){
00093         int nnodes = FEM_Mesh_get_length(meshID,nodeID);
00094         int nelems = FEM_Mesh_get_length(meshID,elemID);
00095 
00096         printf("%d %d \n",nnodes,nelems);       
00097         for(int k=0;k<nnodes;k++){
00098                 printf(" node %d ( %.6f %.6f )\n",k,coord[2*k+0],coord[2*k+1]);
00099         }
00100         REFINE2D_Split(nnodes,coord,nelems,desiredAreas);
00101         int nSplits=REFINE2D_Get_Split_Length();
00102         printf("called REFINE2D_Split nSplits %d\n",nSplits);
00103         
00104   
00105         if(nSplits == 0){
00106                 return;
00107         }
00108 
00109         /*Copy the cordinates of the nodes into a vector, 
00110                 the cordinates of the new nodes will be inserted
00111                 into this vector and will be used to sort all the
00112                 nodes on the basis of the distance from origin
00113         */
00114         CkVec<double> coordVec;
00115         for(int i=0;i<nnodes*2;i++){
00116                 coordVec.push_back(coord[i]);
00117         }
00118         
00119         /*find out the attributes of the node 
00120         */
00121         FEM_Entity *e=FEM_Entity_lookup(meshID,nodeID,"REFINE2D_Mesh");
00122         CkVec<FEM_Attribute *> *attrs = e->getAttrVec();
00123         
00124         /*
00125                 Get the connectivity table of the elements
00126         */
00127         FEM_Entity *elem = FEM_Entity_lookup(meshID,elemID,"REFIN2D_Mesh_elem");
00128         CkVec<FEM_Attribute *> *elemattrs = elem->getAttrVec();
00129         FEM_Attribute *connAttr = elem->lookup(FEM_CONN,"REFINE2D_Mesh");
00130         if(connAttr == NULL){
00131                 CkAbort("Grrrr element without connectivity \n");
00132         }
00133         AllocTable2d<int> &connTable = ((FEM_IndexAttribute *)connAttr)->get();
00134 
00135         /*
00136                 Get the FEM_BOUNDARY data of sparse elements and load it into a hashtable
00137                 indexed by the 2 node ids that make up the edge. The data in the hashtable
00138                 is the index number of the sparse element
00139         */
00140         FEM_Entity *sparse;
00141         CkVec<FEM_Attribute *> *sparseattrs;
00142         FEM_Attribute *sparseConnAttr, *sparseBoundaryAttr;
00143         AllocTable2d<int> *sparseConnTable, *sparseBoundaryTable;
00144         CkHashtableT<intdual,int> nodes2sparse;
00145         if(sparseID != -1){
00146                 sparse = FEM_Entity_lookup(meshID,sparseID,"REFINE2D_Mesh_sparse");
00147                 sparseattrs = sparse->getAttrVec();
00148                 sparseConnAttr = sparse->lookup(FEM_CONN,"REFINE2D_Mesh_sparse");
00149                 sparseConnTable = &(((FEM_IndexAttribute *)sparseConnAttr)->get());
00150                 sparseBoundaryAttr = sparse->lookup(FEM_BOUNDARY,"REFINE2D_Mesh_sparse");
00151                 if(sparseBoundaryAttr == NULL){
00152                          CkAbort("Specified sparse elements without boundary conditions");
00153                 }
00154                 /*
00155                         since the default value in the hashtable is 0, to 
00156                         distinguish between uninserted keys and the sparse element
00157                         with index 0, the index of the sparse elements is incremented
00158                         by 1 while inserting.
00159                 */
00160 //              printf("[%d] Sparse elements\n",FEM_My_partition());
00161                 for(int j=0;j<sparse->size();j++){
00162                         int *cdata = (*sparseConnTable)[j];
00163         //              printf("%d < %d,%d > \n",j,cdata[0],cdata[1]);
00164                         nodes2sparse.put(intdual(cdata[0],cdata[1])) = j+1;
00165                 }
00166         }
00167         
00168   for (int splitNo=0;splitNo<nSplits;splitNo++){
00169     int tri,A,B,C,D;
00170     double frac;
00171                 // current number of nodes in the mesh
00172                 int cur_nodes = FEM_Mesh_get_length(meshID,nodeID);
00173                 int *connData = connTable.getData();
00174                 int flags;
00175 
00176 
00177                 REFINE2D_Get_Split(splitNo,(int *)(connData),&tri,&A,&B,&C,&frac,&flags);
00178                 if((flags & 0x1) || (flags & 0x2)){                     
00179                         //new node 
00180                         D = cur_nodes;
00181       CkPrintf("---- Adding node %d\n",D);                                      
00182                 /*      lastA=A;
00183                         lastB=B;
00184                         lastD=D;*/
00185       if (A>=cur_nodes) CkAbort("Calculated A is invalid!");
00186       if (B>=cur_nodes) CkAbort("Calculated B is invalid!");
00187                         CmiMemoryCheck();
00188                         /*
00189                                 set the data values of the new node
00190                         */
00191                         printf("new node added %d between %d (%.6f,%.6f) %d (%.6f,%.6f)\n",D,A,coordVec[2*A],coordVec[2*A+1],B,coordVec[2*B],coordVec[2*B+1]);                                          
00192                         e->setLength(cur_nodes+1);
00193                         for(int i=0;i<attrs->size();i++){
00194                                 FEM_Attribute *a = (FEM_Attribute *)(*attrs)[i];
00195                                 if(a->getAttr()<FEM_ATTRIB_TAG_MAX){
00196                                         FEM_DataAttribute *d = (FEM_DataAttribute *)a;
00197                                         d->interpolate(A,B,D,frac);
00198                                 }else{
00199                                                 /*The boundary value of a new node should be the 
00200                                                 boundary value of the edge(sparse element) that contains
00201                                                 the two nodes */
00202                                                 if(a->getAttr() == FEM_BOUNDARY){
00203                                                         if(sparseID != -1){
00204                                                                 int sidx = nodes2sparse.get(intdual(A,B))-1;
00205                                                                 if(sidx == -1){
00206                                                                         CkAbort("no sparse element between these 2 nodes, are they really connected ??");
00207                                                                 }
00208                                                                 sparseBoundaryTable = &(((FEM_DataAttribute *)sparseBoundaryAttr)->getInt());
00209                                                                 int boundaryVal = ((*sparseBoundaryTable)[sidx])[0];
00210                                                                 (((FEM_DataAttribute *)a)->getInt()[D])[0] = boundaryVal;
00211                                                         }else{
00212                                                                 /*
00213                                                                         if sparse elements don't exist then just do simple
00214                                                                         interpolation
00215                                                                 */
00216                                                                 FEM_DataAttribute *d = (FEM_DataAttribute *)a;
00217                                                                 d->interpolate(A,B,D,frac);
00218                                                         }
00219                                                 }
00220                                 }       
00221                         }
00222                         int AandB[2];
00223       AandB[0]=A;
00224                   AandB[1]=B;
00225       /* Add a new node D between A and B in the communication list*/
00226                   IDXL_Add_entity(FEM_Comm_shared(meshID,nodeID),D,2,AandB);
00227                         double Dx = coord[2*A]*(1-frac)+frac*coord[2*B];
00228                         double Dy = coord[2*A+1]*(1-frac)+frac*coord[2*B+1];                            
00229                         coordVec.push_back(Dx);
00230                         coordVec.push_back(Dy);
00231                         /*
00232                                 add the new sparse element <D,B> and modify the connectivity of the old one
00233                                 from <A,B> to <A,D> and change the hashtable to reflect that change
00234                         */
00235                         if(sparseID != -1){
00236                                 int oldsidx = nodes2sparse.get(intdual(A,B))-1;
00237                                 int newsidx = sparse->size();
00238                                 sparse->setLength(newsidx+1);
00239                                 for(int satt = 0;satt<sparseattrs->size();satt++){
00240                                         if((*sparseattrs)[satt]->getAttr() == FEM_CONN){
00241                                                 /*
00242                                                         change the conn of the old sparse to A,D
00243                                                         and new one to B,D
00244                                                 */
00245                                                 sparseConnTable = &(((FEM_IndexAttribute *)sparseConnAttr)->get());
00246                                                 int *oldconn = (*sparseConnTable)[oldsidx];
00247                                                 int *newconn = (*sparseConnTable)[newsidx];
00248                                                 oldconn[0] = A;
00249                                                 oldconn[1] = D;
00250                                                 
00251                                                 newconn[0] = D;
00252                                                 newconn[1] = B;
00253                                                 
00254 //                                              printf("<%d,%d> edge being split into <%d,%d> <%d,%d> \n",A,B,A,D,D,B);
00255                                         }else{
00256                                                 /*
00257                                                         apart from conn copy everything else
00258                                                 */
00259                                                 FEM_Attribute *attr = (FEM_Attribute *)(*sparseattrs)[satt];
00260                                                 attr->copyEntity(newsidx,*attr,oldsidx);
00261                                         }
00262                                 }
00263                                 /*
00264                                         modify the hashtable - delete the old edge
00265                                         and the new ones
00266                                 */
00267                                 nodes2sparse.remove(intdual(A,B));
00268                                 nodes2sparse.put(intdual(A,D)) = oldsidx+1;
00269                                 nodes2sparse.put(intdual(D,B)) = newsidx+1;
00270                         }
00271                         
00272                 }
00273                 
00274                 
00275                 
00276                 //add a new triangle
00277                 /*TODO: replace  FEM_ELEM with parameter*/
00278                 int newTri =  FEM_Mesh_get_length(meshID,elemID);
00279     CkPrintf("---- Adding triangle %d after splitting %d \n",newTri,tri);
00280                 elem->setLength(newTri+1);
00281                 for(int j=0;j<elemattrs->size();j++){
00282                         if((*elemattrs)[j]->getAttr() == FEM_CONN){
00283                                 CkPrintf("elem attr conn code %d \n",(*elemattrs)[j]->getAttr());
00284                                 //it is a connectivity attribute.. get the connectivity right
00285                                 FEM_IndexAttribute *connAttr = (FEM_IndexAttribute *)(*elemattrs)[j];
00286                                 AllocTable2d<int> &table = connAttr->get();
00287                                 CkPrintf("Table of connectivity attribute starts at %p width %d \n",table[0],connAttr->getWidth());
00288                                 int *oldRow = table[tri];
00289                                 int *newRow = table[newTri];
00290                                 for (int i=0;i<3;i++){
00291                       if (oldRow[i]==A){
00292                                                 oldRow[i]=D;    
00293                                                 CkPrintf("In triangle %d %d replaced by %d \n",tri,A,D);
00294                                         }       
00295                                 }       
00296                                 for (int i=0; i<3; i++) {
00297                       if (oldRow[i] == B){
00298                                                 newRow[i] = D;
00299                                         }       
00300                 else if (oldRow[i] == C){
00301                                                 newRow[i] = C;
00302                                         }       
00303                       else if (oldRow[i] == D){
00304                                                 newRow[i] = A;
00305                                         }       
00306                         }
00307                                 CkPrintf("New Triangle %d  (%d %d %d) conn %p\n",newTri,newRow[0],newRow[1],newRow[2],newRow);
00308                         }else{
00309                                 FEM_Attribute *elattr = (FEM_Attribute *)(*elemattrs)[j];
00310                                 if(elattr->getAttr() < FEM_ATTRIB_FIRST){ 
00311                                         elattr->copyEntity(newTri,*elattr,tri);
00312                                 }       
00313                         }
00314                 }
00315                 if(sparseID != -1){
00316                         /*
00317                                 add the sparse element (edge between C and D)
00318                         */
00319                         int cdidx = sparse->size();
00320                         sparse->setLength(cdidx+1);
00321                         for(int satt = 0; satt < sparseattrs->size();satt++){
00322                                         if((*sparseattrs)[satt]->getAttr() == FEM_CONN){
00323                                                 sparseConnTable = &(((FEM_IndexAttribute *)sparseConnAttr)->get());
00324                                                 int *cdconn = (*sparseConnTable)[cdidx];
00325                                                 cdconn[0]=C;
00326                                                 cdconn[1]=D;
00327                                         }
00328                                         if((*sparseattrs)[satt]->getAttr() == FEM_BOUNDARY){
00329                                                 /*
00330                                                         An edge connecting C and D has to be an internal edge
00331                                                 */
00332                                                 sparseBoundaryTable = &(((FEM_DataAttribute *)sparseBoundaryAttr)->getInt());
00333                                                 ((*sparseBoundaryTable)[cdidx])[0] = 0;
00334                                         }
00335                         }
00336                         nodes2sparse.put(intdual(C,D)) = cdidx+1;
00337                 }
00338                 
00339         }
00340         printf("Cordinate list length %d \n",coordVec.size()/2);
00341         IDXL_Sort_2d(FEM_Comm_shared(meshID,nodeID),coordVec.getVec());
00342         
00343 }
00344 
00345 FDECL void FTN_NAME(FEM_REFINE2D_SPLIT,fem_refine2d_split)(int *meshID,int *nodeID,double *coord,int *elemID,double *desiredAreas){
00346         FEM_REFINE2D_Split(*meshID,*nodeID,coord,*elemID,desiredAreas);
00347 }
libs/ck-libs/tmr/femrefine.C
