Charm++: Charm Source Code Documentation

00001 // Simple sample implementation of refinement interface.
00002 // Orion Sky Lawlor, olawlor@acm.org, 4/9/2002
00003 // Modified by Terry Wilmarth, wilmarth@cse.uiuc.edu, 4/16/2002
00004 #include <stdlib.h>
00005 #include <unistd.h>
00006 #include <stdio.h>
00007 #include <math.h>
00008 #include "charm++.h"
00009 #include "charm-api.h"
00010 #include "tcharm.h"
00011 #include "mpi.h"
00012 #include "tri.h"
00013 #include "refine.h"
00014 
00015 /********************* Attach *****************/
00016 CLINKAGE void REFINE2D_Init(void) {
00017   TCHARM_API_TRACE("REFINE2D_Init", "refine");
00018   TCharm *tc=TCharm::get();
00019   
00020   // Create the refinement array
00021   MPI_Comm comm=MPI_COMM_WORLD; /* fixme: make user pass in communicator */
00022   int rank; MPI_Comm_rank(comm,&rank);
00023   CkArrayID refArrayID;
00024   if (rank==0) 
00025   { /* Master creates the (empty) array and broadcasts arrayID: */
00026     CkArrayOptions opts;
00027     opts.bindTo(tc->getProxy());
00028     refArrayID=CProxy_chunk::ckNew(new chunkMsg, opts);
00029   }
00030   MPI_Bcast(&refArrayID,sizeof(refArrayID),MPI_BYTE, 0,comm);
00031   mesh=refArrayID; /* set up readonly proxy */
00032   
00033   // Now everybody inserts their element into the array
00034   chunkMsg *cm = new chunkMsg;
00035   cm->nChunks = tc->getNumElements();
00036   cm->myThreads = tc->getProxy();
00037   mesh[rank].insert(cm);
00038   tc->suspend(); /* will resume from chunk constructor */
00039 }
00040 FLINKAGE void FTN_NAME(REFINE2D_INIT,refine2d_init)(void)
00041 {
00042   REFINE2D_Init();
00043 }
00044 
00045 /******************** NewMesh *******************/
00046 CLINKAGE void REFINE2D_NewMesh(int meshID,int nEl,int nGhost,int nnodes,const int *conn,const int *gid,const int *boundaries, const int *edgeBounds, const int *edgeConn, int nEdges)
00047 {
00048   TCHARM_API_TRACE("REFINE2D_NewMesh", "refine");
00049   if (!CtvAccess(_refineChunk))
00050     CkAbort("Forgot to call REFINE_Attach!!\n");
00051     
00052   CtvAccess(_refineChunk)->newMesh(meshID, nEl, nGhost, conn, gid, nnodes, 
00053                    boundaries, nEdges, edgeConn, edgeBounds, 0);
00054   MPI_Barrier(MPI_COMM_WORLD);
00055   CkWaitQD();
00056 }
00057 /*
00058 FLINKAGE void FTN_NAME(REFINE2D_NEWMESH,refine2d_newmesh)
00059 (int *nEl,int *nGhost,int nnodes,const int *conn,const int *gid,const int *boundaries, const int **edgeBoundaries)
00060 {
00061   TCHARM_API_TRACE("REFINE2D_NewMesh", "refine");
00062   if (!CtvAccess(_refineChunk))
00063     CkAbort("Forgot to call REFINE_Attach!!\n"); 
00064   
00065   CtvAccess(_refineChunk)->newMesh(*nEl, *nGhost,conn, gid, nnodes, 
00066                                    boundaries, edgeBoundaries, 1);
00067   MPI_Barrier(MPI_COMM_WORLD);
00068   CkWaitQD();
00069 }*/
00070 
00071 /********************** Splitting ******************/
00072 class refineResults {
00073   std::vector<refineData> res;
00074   int nResults;
00075   
00076   //Return anything on [0,2] than a or b
00077   int otherThan(int a,int b) {
00078     if (a==b) CkAbort("Opposite node is moving!");
00079     for (int i=0;i<3;i++)
00080       if (i!=a && i!=b) return i;
00081     CkAbort("Logic error in refine.C::otherThan");
00082     return -1;
00083   }
00084 public:
00085   refineResults(void) {nResults=0;}
00086     refineData createRefineData(int tri, int A, int B, int C,  
00087         int D, int _new, double frac,int flag, int origEdgeB, int newEdge1B, 
00088         int newEdge2B){
00089         refineData d;
00090         d.tri = tri;
00091         d.A = A;
00092         d.B = B;
00093         d.C = C;
00094         d.D = D;
00095         d._new = _new;
00096         d.frac = frac;
00097         d.flag = flag;
00098         d.origEdgeB = origEdgeB;
00099         d.newEdge1B = newEdge1B;
00100         d.newEdge2B = newEdge2B;
00101         return d;
00102     }
00103     void add(refineData &d){
00104         nResults++;
00105         res.push_back(d);
00106     };
00107     
00108   int countResults(void) const {return nResults;}
00109     
00110   void extract(int i, refineData *d) {
00111         *d =  res[i];   
00112   }
00113 };
00114 void FEM_Modify_IDXL(FEM_Refine_Operation_Data *data,refineData &d);
00115 
00116 class resultsRefineClient : public refineClient {
00117   refineResults *res;
00118     FEM_Refine_Operation_Data *data;
00119 public:
00120   resultsRefineClient(refineResults *res_,FEM_Refine_Operation_Data *data_) :res(res_),data(data_) {}
00121  /* void split(int tri, int side, int node, double frac) {
00122 #if 0
00123     //Convert from "tri.C edges" to sensible edges
00124     if (side==1) side=2;
00125     else if (side==2) side=1;
00126 #endif
00127     res->add(tri, side, node, frac);
00128   }*/
00129     //for the explanation of A,B,C,D look at diagram in refine.h
00130   void split(int tri, int A, int B, int C,  
00131          int D, int _new, double frac,int flag, int origEdgeB, int newEdge1B, 
00132          int newEdge2B) {
00133     refineData d = res->createRefineData(tri,A,B,C,D, _new,frac,flag,origEdgeB,
00134         newEdge1B,newEdge2B);
00135         
00136         FEM_Modify_IDXL(data,d);
00137 
00138         res->add(d);
00139   }
00140 };
00141 
00142 class coarsenResults {
00143   // coarsenData is defined in refine.h
00144   std::vector<coarsenData> res;
00145 public:
00146   coarsenResults(){}
00147   coarsenData addCollapse(int elementID, int nodeToKeep, int nodeToDelete,
00148               double nX, double nY, int flag, int boundFlag, 
00149               double frac)
00150   {
00151     coarsenData d;
00152     d.type = COLLAPSE;
00153     d.data.cdata.elemID = elementID;
00154     d.data.cdata.nodeToKeep = nodeToKeep;
00155     d.data.cdata.nodeToDelete = nodeToDelete;
00156     d.data.cdata.newX = nX;
00157     d.data.cdata.newY = nY;
00158     d.data.cdata.flag = flag;
00159     d.data.cdata.boundaryFlag = boundFlag;
00160     d.data.cdata.frac = frac;
00161     return d;
00162   };
00163   coarsenData addUpdate(int nodeID, double newX, double newY, int boundaryFlag)
00164   {
00165     coarsenData d;
00166     d.type = UPDATE;
00167     d.data.udata.nodeID = nodeID;
00168     d.data.udata.newX = newX;
00169     d.data.udata.newY = newY;
00170     d.data.udata.boundaryFlag = boundaryFlag;
00171     //      res.push_back(d);
00172     return d;
00173   };
00174   coarsenData addReplaceDelete(int elemID, int relnodeID, int oldNodeID,
00175                    int newNodeID)
00176   {
00177     coarsenData d;
00178     d.type = REPLACE;
00179     d.data.rddata.elemID = elemID;
00180     d.data.rddata.relnodeID = relnodeID;
00181     d.data.rddata.oldNodeID = oldNodeID;
00182     d.data.rddata.newNodeID = newNodeID;
00183     //      res.push_back(d);
00184     return d;
00185   };
00186   int countResults(){return res.size();}
00187   /*void extract(int i,int *conn,int *tri,int *nodeToThrow,int *nodeToKeep,double *nx,double *ny,int *flag,int idxbase){
00188     int t;
00189     t = res[i].elemID;
00190     *tri = t +idxbase;
00191     *nodeToKeep = conn[3*t+res[i].nodeToKeep];
00192     int n1 = res[i].collapseEdge;
00193     int n2 = (res[i].collapseEdge+1)%3;
00194     if(res[i].nodeToKeep == n1){
00195     *nodeToThrow = conn[3*t+n2];
00196     }else{
00197     *nodeToThrow = conn[3*t+n1];
00198     }
00199     *nx = res[i].nx;
00200     *ny = res[i].ny;
00201     *flag = res[i].flag;
00202     }*/
00203   void extract(int i, coarsenData *output){
00204     *output = res[i];
00205   }
00206 };
00207 
00208 // Modifying the code so that instead of being stored the results are
00209 // processed immediately. The coarsen client is used simply to create
00210 // and return an appropriate coarsenData structure.
00211 class FEM_Operation_Data;
00212 void FEM_Coarsen_Operation(FEM_Operation_Data *coarsen_data, 
00213                coarsenData &operation);
00214 
00215 class resultsCoarsenClient : public refineClient {
00216   coarsenResults *res;
00217   FEM_Operation_Data *data;
00218 public:
00219   resultsCoarsenClient(coarsenResults *res_, FEM_Operation_Data *data_=NULL) 
00220     : res(res_),data(data_){};
00221   void collapse(int elementID, int nodeToKeep, int nodeToDelete, double nX,
00222         double nY, int flag, int b, double frac)
00223   {
00224     coarsenData d = res->addCollapse(elementID, nodeToKeep, nodeToDelete, nX,
00225                      nY, flag, b, frac);
00226     FEM_Coarsen_Operation(data,d);
00227   }
00228   void nodeUpdate(int nodeID, double newX, double newY, int boundaryFlag)
00229   {
00230     coarsenData d = res->addUpdate(nodeID,newX,newY,boundaryFlag);
00231     FEM_Coarsen_Operation(data,d);
00232   }
00233   void nodeReplaceDelete(int elementID, int relnodeID, int oldNodeID, 
00234              int newNodeID)
00235   {
00236     coarsenData d = res->addReplaceDelete(elementID, relnodeID, oldNodeID,
00237                       newNodeID);
00238     FEM_Coarsen_Operation(data,d);
00239   }
00240 };
00241 
00242 
00243 // this function should be called from a thread
00244 CLINKAGE void REFINE2D_Split(int nNode,double *coord,int nEl,double *desiredArea,FEM_Refine_Operation_Data *refine_data)
00245 {
00246   TCHARM_API_TRACE("REFINE2D_Split", "refine");
00247   chunk *C = CtvAccess(_refineChunk);
00248   if (!C)
00249     CkAbort("REFINE2D_Split failed> Did you forget to call REFINE2D_Attach?");
00250   C->refineResultsStorage=new refineResults;
00251   resultsRefineClient client(C->refineResultsStorage,refine_data);
00252 
00253   C->updateNodeCoords(nNode, coord, nEl);
00254   C->multipleRefine(desiredArea, &client);
00255 }
00256 
00257 CLINKAGE void REFINE2D_Coarsen(int nNode, double *coord, int nEl,
00258                 double *desiredArea, FEM_Operation_Data *data)
00259 {
00260   TCHARM_API_TRACE("REFINE2D_Coarsen", "coarsen");
00261   chunk *C = CtvAccess(_refineChunk);
00262   if (!C)
00263     CkAbort("REFINE2D_Split failed> Did you forget to call REFINE2D_Attach?");
00264   C->coarsenResultsStorage=new coarsenResults;
00265   resultsCoarsenClient client(C->coarsenResultsStorage,data);
00266 
00267   C->updateNodeCoords(nNode, coord, nEl);
00268   C->multipleCoarsen(desiredArea, &client);
00269 }
00270 
00271 
00272 FLINKAGE void FTN_NAME(REFINE2D_SPLIT,refine2d_split)
00273    (int *nNode,double *coord,int *nEl,double *desiredArea,FEM_Refine_Operation_Data *data)
00274 {
00275   REFINE2D_Split(*nNode,coord,*nEl,desiredArea,data);
00276 }
00277 
00278 static refineResults *getResults(void) {
00279   chunk *C = CtvAccess(_refineChunk);
00280   if (!C)
00281     CkAbort("Did you forget to call REFINE2D_Init?");
00282   refineResults *ret=C->refineResultsStorage;
00283   if (ret==NULL)
00284     CkAbort("Did you forget to call REFINE2D_Begin?");
00285   return ret;
00286 }
00287 
00288 CLINKAGE int REFINE2D_Get_Split_Length(void)
00289 {
00290   TCHARM_API_TRACE("REFINE2D_Get_Split_Length", "refine");
00291   return getResults()->countResults();
00292 }
00293 FLINKAGE int FTN_NAME(REFINE2D_GET_SPLIT_LENGTH,refine2d_get_split_length)(void)
00294 {
00295   return REFINE2D_Get_Split_Length();
00296 }
00297 
00298 CLINKAGE void REFINE2D_Get_Split
00299     (int splitNo,refineData *d)
00300 {
00301   TCHARM_API_TRACE("REFINE2D_Get_Split", "refine");
00302   refineResults *r=getResults();
00303   r->extract(splitNo,d);
00304 }
00305 FLINKAGE void FTN_NAME(REFINE2D_GET_SPLIT,refine2d_get_split)
00306     (int *splitNo,refineData *d)
00307 {
00308   TCHARM_API_TRACE("REFINE2D_Get_Split", "refine");
00309   refineResults *r=getResults();
00310   r->extract(*splitNo-1,d);
00311 }
00312 
00313 static coarsenResults *getCoarsenResults(void) {
00314   chunk *C = CtvAccess(_refineChunk);
00315   if (!C)
00316     CkAbort("Did you forget to call REFINE2D_Init?");
00317   coarsenResults *ret=C->coarsenResultsStorage;
00318   if (ret==NULL)
00319     CkAbort("Did you forget to call REFINE2D_Coarsen?");
00320   return ret;
00321 }
00322 
00323 
00324 CLINKAGE int REFINE2D_Get_Collapse_Length(){
00325   return getCoarsenResults()->countResults();
00326 }
00327 /*
00328 CLINKAGE void REFINE2D_Get_Collapse(int i,int *conn,int *tri,int *nodeToThrow,int *nodeToKeep,double *nx,double *ny,int *flag,int idxbase){
00329     return getCoarsenResults()->extract(i,conn,tri,nodeToThrow,nodeToKeep,nx,ny,flag,idxbase);
00330 }*/
00331 
00332 CLINKAGE void REFINE2D_Get_Collapse(int i,coarsenData *output){
00333     getCoarsenResults()->extract(i,output);
00334 }
00335 
00336 /********************* Check *****************/
00337 static int checkElement(chunk *C,const element &e,const int *uc,int idxBase)
00338 {
00339   int nMismatch=0;
00340   //Make sure the nodes are equal:
00341   /* This check no longer applies...
00342   for (int j=0;j<3;j++) {
00343     if (e.nodes[j].cid!=C->cid)
00344         CkAbort("REFINE2D_Check> Triangle has non-local node!\n");
00345     if (uc[j]-idxBase!=e.nodes[j].idx) {
00346       //Very bad-- mismatching indices
00347       nMismatch++;
00348     }
00349   }
00350   */
00351   return nMismatch;
00352 }
00353 
00354 static void checkConn(int nEl,const int *conn,int idxBase,int nNode)
00355 {
00356   chunk *C = CtvAccess(_refineChunk);
00357   if (!C) CkAbort("Did you forget to call REFINE2D_Attach?");
00358   //This check stolen from updateNodeCoords:
00359   if (nEl != C->numElements || nNode != C->numNodes) {
00360     CkPrintf("ERROR: inconsistency in REFINE2D_Check on chunk %d:\n"
00361        "  your nEl (%d); my numElements (%d)\n"
00362        "  your nNode (%d); my numNodes (%d)\n",
00363        C->cid, nEl, C->numElements, nNode,C->numNodes);
00364     CkAbort("User code/library numbering inconsistency in REFINE2D");
00365   }
00366   int i;
00367   int nErrs=0;
00368   for (i=0;i<nEl;i++) {
00369     const element &e=C->theElements[i];
00370     const int *uc=&conn[3*i];
00371     int elErrs=checkElement(C,e,uc,idxBase);
00372     nErrs+=elErrs;
00373     if (elErrs!=0) { //A bad element-- print out debugging info:
00374       /* This check is no longer valid...
00375      CkError("REFINE2D Chunk %d>   Your triangle %d: %d %d %d\n"
00376      "REFINE2D Chunk %d> but my triangle %d: %d %d %d\n",
00377      C->cid+idxBase,i+idxBase, uc[0], uc[1], uc[2],
00378      C->cid+idxBase,i+idxBase,
00379      e.nodes[0].idx+idxBase,e.nodes[1].idx+idxBase,e.nodes[2].idx+idxBase);
00380       */
00381     }
00382   }
00383   if (nErrs!=0) {
00384     CkAbort("REFINE2D_Check> Major errors found. Exiting.");
00385   }
00386 }
00387 
00388 CLINKAGE void REFINE2D_Check(int nEl,const int *conn,int nNodes) {
00389   checkConn(nEl,conn,0,nNodes);
00390 }
00391 FLINKAGE void FTN_NAME(REFINE2D_CHECK,refine2d_check)
00392   (int *nEl,const int *conn,int *nNodes)
00393 {
00394   checkConn(*nEl,conn,1,*nNodes);
00395 }
00396 
00397 
00398
libs/ck-libs/TMRC2D/refine.C
