Charm++: Charm Source Code Documentation

00001 // Triangular Mesh Refinement Framework - 2D (TMR)
00002 // Created by: Terry L. Wilmarth
00003 
00004 #include <stdlib.h>
00005 #include <stdio.h>
00006 #include "tri.h"
00007 
00008 // I will keep this code until I die!
00009 //#define accessLock() {printf("%s : %d calling accessLock \n",__FILE__,__LINE__); accessLock1();}
00010 //readonlys
00011 CProxy_chunk mesh;
00012 CtvDeclare(chunk *, _refineChunk);
00013 
00014 void refineChunkInit(void) {
00015   CtvInitialize(chunk *, _refineChunk);
00016 }
00017 
00018 chunk::chunk(chunkMsg *m)
00019   : TCharmClient1D(m->myThreads), sizeElements(0), sizeEdges(0), sizeNodes(0),
00020     firstFreeElement(0), firstFreeEdge(0), firstFreeNode(0),
00021     edgesSent(0), edgesRecvd(0), first(0),
00022     coarsenElements(NULL), refineElements(NULL), refineStack(NULL),
00023     refineHeapSize(0), coarsenHeapSize(0), refineTop(0),
00024     additions(0), debug_counter(0), refineInProgress(0), coarsenInProgress(0),
00025     meshLock(0), meshExpandFlag(0), 
00026     numElements(0), numEdges(0), numNodes(0), numGhosts(0), theClient(NULL),
00027     elementSlots(0), edgeSlots(0), nodeSlots(0), lock(0), lockCount(0),
00028     lockHolderIdx(-1), lockHolderCid(-1), lockPrio(-1.0), lockList(NULL)
00029 {
00030   refineResultsStorage=NULL;
00031   cid = thisIndex;
00032   numChunks = m->nChunks;
00033   CkFreeMsg(m);
00034   tcharmClientInit();
00035   thread->resume();
00036 }
00037 
00038 void chunk::addRemoteEdge(int elem, int localEdge, edgeRef er)
00039 {
00040   accessLock();
00041   CkAssert(localEdge >=0);
00042   CkAssert(localEdge < 3);
00043   CkAssert(er.cid > -1);
00044   CkAssert(er.idx > -1);
00045   if ((theElements[elem].edges[localEdge].cid > -1) &&
00046       (theElements[elem].edges[localEdge].idx > -1)){
00047     CkPrintf("TMRC2D: [%d] WARNING: addRemoteEdge replacing non-null edgeRef!\n", cid);
00048   } 
00049   theElements[elem].set(localEdge, er);
00050 #ifdef TDEBUG3
00051   CkPrintf("TMRC2D: [%d] addRemoteEdge on element %d", cid, elem);
00052 #endif
00053   edgesRecvd++;
00054   releaseLock();
00055 }
00056 
00057 void chunk::refineElement(int idx, double area)
00058 { // Reduce element's targetArea to indicate need for refinement
00059   if (!theElements[idx].isPresent()) return;
00060   accessLock();
00061   theElements[idx].resetTargetArea(area);
00062   refineStack[refineTop].elID = idx;
00063   refineStack[refineTop].len = -1.0;
00064   refineTop++;
00065   int pos = refineHeapSize;
00066   while (pos >= 1) {
00067     if (refineElements[pos].elID == idx)
00068       refineElements[pos].elID = -1;
00069     pos--;
00070   }
00071   releaseLock();
00072   if (!refineInProgress) { // if refine loop not running
00073     refineInProgress = 1;
00074     mesh[cid].refiningElements(); // start it up
00075   }
00076 }
00077 
00078 void chunk::refiningElements()
00079 {
00080   int i;
00081   while (refineHeapSize>0 || refineTop > 0) { // loop through the elements
00082     if (refineTop>0) {
00083       refineTop--;
00084       i=refineStack[refineTop].elID;
00085     }
00086     else
00087       i=Delete_Min(0);
00088     if ((i != -1) && theElements[i].isPresent()) {
00089       // element i has a lower target area -- needs to refine
00090       theElements[i].refine(); // refine the element
00091       adjustMesh();
00092     }
00093     CthYield(); // give other chunks on the same PE a chance
00094   }
00095   sanityCheck(); // quietly make sure mesh is in shape
00096   refineInProgress = 0; // nothing needs refinement; turn refine loop off
00097   //    if (CkMyPe() == 0) for (int j=0; j<5; j++) mesh[j].print();
00098 }
00099 
00100 
00101 void chunk::coarsenElement(int idx, double area)
00102 { 
00103   if (!theElements[idx].isPresent()) return;
00104   double theArea;
00105   accessLock();
00106   theArea = theElements[idx].getArea();
00107   if (area > theElements[idx].getTargetArea())
00108     theElements[idx].resetTargetArea(area);
00109   Insert(idx, -1.0, 1);
00110   int pos = coarsenHeapSize;
00111   while (pos > 1) {
00112     if (coarsenElements[pos].elID == idx)
00113       coarsenElements[pos].elID = -1;
00114     pos--;
00115   }
00116   releaseLock();
00117   if (!coarsenInProgress) { // if coarsen loop not running
00118     coarsenInProgress = 1;
00119     mesh[cid].coarseningElements(); // start it up
00120   }
00121 }
00122 
00123 void chunk::coarseningElements()
00124 {
00125   int i;
00126   double area, qFactor, targetArea;
00127   while (coarsenHeapSize > 0) { // loop through the elements
00128     i=Delete_Min(1);
00129     CkPrintf("Loop at element %d\n", i);
00130     if ((i != -1) && theElements[i].isPresent()) {
00131       CkPrintf("Checking element %d\n", i);
00132       area = theElements[i].getArea();
00133       targetArea = theElements[i].getTargetArea();
00134       qFactor = theElements[i].getAreaQuality();
00135       if ((theElements[i].getTargetArea()*COARSEN_PRECISION > area) || 
00136       (area == 0.0) || (qFactor < QUALITY_MIN)) { 
00137     CkPrintf("Coarsening element %d\n", i);
00138     theElements[i].coarsen(); // coarsen the element
00139       }
00140     }
00141     CthYield(); // give other chunks on the same PE a chance
00142   }
00143   coarsenInProgress = 0;  // turn coarsen loop off
00144   if (coarsenElements) delete [] coarsenElements;
00145   coarsenElements = new elemHeap[numElements+1];
00146   coarsenElements[0].elID=-1;
00147   coarsenElements[0].len=-2.0;
00148   for (i=0; i<elementSlots; i++) {
00149     if (theElements[i].isPresent() && (theElements[i].nonCoarsenCount<1)) {
00150       area = theElements[i].getArea();
00151       if (area == 0.0) {
00152     CkPrintf("Element[%d] has area %1.10e!\n", i, area);
00153       }
00154       targetArea = theElements[i].getTargetArea();
00155       double shortEdgeLen, angle;
00156       shortEdgeLen=theElements[i].getShortestEdge(&angle);
00157       qFactor = theElements[i].getAreaQuality();
00158       if ((targetArea*COARSEN_PRECISION > area) || (qFactor < QUALITY_MIN)) {
00159     CkPrintf("Element[%d] has area %1.10e target %1.10e qFactor %1.10e\n", i, area, targetArea, qFactor);
00160     Insert(i, shortEdgeLen*qFactor, 1);
00161       }
00162     }
00163   }
00164   sanityCheck(); // quietly make sure mesh is in shape
00165   if ((coarsenHeapSize>0) && !coarsenInProgress) {
00166     coarsenInProgress = 1;
00167     mesh[cid].coarseningElements();
00168   }
00169 }
00170 
00171 // many remote access methods follow
00172 intMsg *chunk::safeToMoveNode(int idx, double x, double y)
00173 {
00174   node foo(x, y);
00175   intMsg *im = new intMsg;
00176   accessLock();
00177   im->anInt = theNodes[idx].safeToMove(foo);
00178   releaseLock();
00179   return im;
00180 }
00181 
00182 splitOutMsg *chunk::split(int idx, elemRef e, int oIdx, int fIdx)
00183 {
00184   splitOutMsg *som = new splitOutMsg;
00185   accessLock();
00186   som->result = theEdges[idx].split(&(som->n), &(som->e), oIdx, fIdx, e, 
00187                     &(som->local), &(som->first), 
00188                     &(som->nullNbr));
00189   releaseLock();
00190   return som;
00191 }
00192 
00193 void chunk::collapse(int idx, elemRef e, int kIdx, int dIdx, elemRef kNbr, 
00194              elemRef dNbr, edgeRef kEdge, edgeRef dEdge, node newN, 
00195              double frac)
00196 {
00197   accessLock();
00198   theEdges[idx].collapse(e, kIdx, dIdx, kNbr, dNbr, kEdge, dEdge, newN, frac);
00199   releaseLock();
00200 }
00201 
00202 splitOutMsg *chunk::flipPreventE(int idx, elemRef e, int kIdx, int dIdx,
00203                  elemRef kNbr, elemRef dNbr, edgeRef kEdge, 
00204                  edgeRef dEdge, node newN)
00205 {
00206   splitOutMsg *som = new splitOutMsg;
00207   accessLock();
00208   som->result = theEdges[idx].flipPrevent(e, kIdx, dIdx, kNbr, dNbr, kEdge, 
00209                       dEdge, newN);
00210   releaseLock();
00211   return som;
00212 }
00213 
00214 void chunk::nodeReplaceDelete(int kIdx, int dIdx, node nn, int shared, 
00215                   int *chk, int *idx)
00216 {
00217   int *jChk, *jIdx;
00218   int jShared;
00219   accessLock();
00220   if (dIdx == -1) { 
00221     if (kIdx != -1) {
00222       theNodes[kIdx].set(nn.X(), nn.Y());
00223       theNodes[kIdx].boundary = nn.boundary;
00224       theNodes[kIdx].fixed = nn.fixed;
00225       jShared = joinCommLists(kIdx, shared, chk, idx, jChk, jIdx);
00226       theClient->nodeUpdate(kIdx, nn.X(), nn.Y(), nn.boundary, jShared, jChk, 
00227                 jIdx);
00228 #ifdef TDEBUG1
00229       CkPrintf("TMRC2D: [%d] (a)theClient->nodeUpdate(%d, %2.10f, %2.10f, %d)\n", cid, kIdx, nn.X(), nn.Y(), nn.boundary);
00230 #endif
00231     }
00232     return;
00233   }
00234   else if (kIdx == -1) {
00235     theNodes[dIdx].set(nn.X(), nn.Y());
00236     theNodes[dIdx].boundary = nn.boundary;
00237     theNodes[dIdx].fixed = nn.fixed;
00238     jShared = joinCommLists(dIdx, shared, chk, idx, jChk, jIdx);
00239     theClient->nodeUpdate(dIdx, nn.X(), nn.Y(), nn.boundary, jShared, jChk, 
00240               jIdx);
00241 #ifdef TDEBUG1
00242     CkPrintf("TMRC2D: [%d] (b)theClient->nodeUpdate(%d, %2.10f, %2.10f, %d)\n", cid, dIdx, nn.X(), nn.Y(), nn.boundary);
00243 #endif
00244   }
00245   else {
00246     removeNode(dIdx);
00247     theNodes[kIdx].set(nn.X(), nn.Y());
00248     theNodes[kIdx].boundary = nn.boundary;
00249     theNodes[kIdx].fixed = nn.fixed;
00250     jShared = joinCommLists(kIdx, shared, chk, idx, jChk, jIdx);
00251     theClient->nodeUpdate(kIdx, nn.X(), nn.Y(), nn.boundary, jShared, jChk, 
00252               jIdx);
00253 #ifdef TDEBUG1
00254     CkPrintf("TMRC2D: [%d] (c)theClient->nodeUpdate(%d, %2.10f, %2.10f, %d)\n", cid, kIdx, nn.X(), nn.Y(), nn.boundary);
00255 #endif
00256     for (int j=0; j<elementSlots; j++) {
00257       if (theElements[j].isPresent()) {
00258     CkAssert(theElements[j].nodes[0] != theElements[j].nodes[1]);
00259     CkAssert(theElements[j].nodes[0] != theElements[j].nodes[2]);
00260     CkAssert(theElements[j].nodes[2] != theElements[j].nodes[1]);
00261     for (int k=0; k<3; k++) {
00262       if (theElements[j].nodes[k] == dIdx) {
00263 #ifdef FLIPTEST
00264         if(theElements[j].nodes[(k+1)%3] == kIdx || theElements[j].nodes[(k+2)%3] == kIdx) break; //don't worry abt this element, it will be deleted
00265         if(theElements[j].flipInverseTest(&(theNodes[dIdx]),&nn)) {
00266           //CkPrintf("Chunk %d, Elem %d is flipping!!\n",cid,j);
00267         }
00268 #endif
00269         theElements[j].nonCoarsenCount = 0;
00270         theElements[j].nodes[k] = kIdx;
00271         theClient->nodeReplaceDelete(j, k, dIdx, kIdx);
00272 #ifdef TDEBUG1
00273         CkPrintf("TMRC2D: [%d] theClient->nodeReplaceDelete(%d, %d, %d, %d)\n", cid, j, k, dIdx, kIdx);
00274 #endif
00275       }
00276       if (theElements[j].nodes[k] == kIdx) {
00277         theElements[j].nonCoarsenCount = 0;
00278       }
00279     }
00280     CkAssert(theElements[j].nodes[0] != theElements[j].nodes[1]);
00281     CkAssert(theElements[j].nodes[0] != theElements[j].nodes[2]);
00282     CkAssert(theElements[j].nodes[2] != theElements[j].nodes[1]);
00283       }
00284     }
00285     for (int j=0; j<edgeSlots; j++) {
00286       if (theEdges[j].isPresent()) {
00287     for (int k=0; k<2; k++) {
00288       if (theEdges[j].nodes[k] == dIdx) {
00289         theEdges[j].nodes[k] = kIdx;
00290 #ifdef TDEBUG1
00291         CkPrintf("TMRC2D: [%d] Edge %d node %d replaced with %d\n", cid, j,
00292              dIdx, kIdx);
00293 #endif
00294       }
00295     }
00296       }
00297     }
00298   }
00299   releaseLock();
00300 }
00301 
00302 boolMsg *chunk::flipPrevent(int kIdx, int dIdx, node nn, int shared, int *chk, int *idx)
00303 {
00304   boolMsg *bm = new boolMsg;
00305   bm->aBool = false;
00306   for (int j=0; j<elementSlots; j++) {
00307     if (theElements[j].isPresent()) {
00308       CkAssert(theElements[j].nodes[0] != theElements[j].nodes[1]);
00309       CkAssert(theElements[j].nodes[1] != theElements[j].nodes[2]);
00310       CkAssert(theElements[j].nodes[2] != theElements[j].nodes[0]);
00311       for (int k=0; k<3; k++) {
00312     if (theElements[j].nodes[k] == dIdx) {
00313       if(theElements[j].nodes[(k+1)%3] == kIdx || theElements[j].nodes[(k+2)%3] == kIdx) break; //don't worry abt this element, it will be deleted
00314       if(theElements[j].flipInverseTest(&(theNodes[dIdx]),&nn)) {
00315         //CkPrintf("Chunk %d, Elem %d is flipping!!\n",cid,j);
00316         bm->aBool = true;
00317         //mesh[*chk].incnonCoarsen(*idx);
00318         return bm;
00319       }
00320     }
00321     if (theElements[j].nodes[k] == kIdx) {
00322       if(theElements[j].nodes[(k+1)%3] == dIdx || theElements[j].nodes[(k+2)%3] == dIdx) break; //don't worry abt this element, it will be deleted
00323       if(theElements[j].flipInverseTest(&(theNodes[kIdx]),&nn)) {
00324         //CkPrintf("Chunk %d, Elem %d is flipping!!\n",cid,j);
00325         bm->aBool = true;
00326         //mesh[*chk].incnonCoarsen(*idx);
00327         return bm;
00328       }
00329     }
00330       }
00331     }
00332   }
00333   return bm;
00334 }
00335 
00336 void chunk::incnonCoarsen(int idx) {
00337   accessLock();
00338   theElements[idx].incnonCoarsen();
00339   releaseLock();
00340   return;
00341 }
00342 
00343 void chunk::resetnonCoarsen(int idx) {
00344   accessLock();
00345   theElements[idx].resetnonCoarsen();
00346   releaseLock();
00347   return;
00348 }
00349 
00350 intMsg *chunk::isPending(int idx, objRef e)
00351 {
00352   intMsg *im = new intMsg;
00353   elemRef eR(e.cid, e.idx);
00354   accessLock();
00355   im->anInt = theEdges[idx].isPending(eR);
00356   releaseLock();
00357   return im;
00358 }
00359 
00360 void chunk::checkPending(int idx, objRef aRef)
00361 {
00362   elemRef eRef;
00363   eRef.idx = aRef.idx; eRef.cid = aRef.cid;
00364   accessLock();
00365   theEdges[idx].checkPending(eRef);
00366   releaseLock();
00367 }
00368 
00369 void chunk::checkPending(int idx, objRef aRef1, objRef aRef2)
00370 {
00371   elemRef eRef1, eRef2;
00372   eRef1.idx = aRef1.idx; eRef1.cid = aRef1.cid;
00373   eRef2.idx = aRef2.idx; eRef2.cid = aRef2.cid;
00374   accessLock();
00375   theEdges[idx].checkPending(eRef1, eRef2);
00376   releaseLock();
00377 }
00378 
00379 void chunk::updateElement(int idx, objRef oldval, objRef newval, int b)
00380 {
00381   elemRef ov, nv;
00382   ov.idx = oldval.idx;   ov.cid = oldval.cid; 
00383   nv.idx = newval.idx;   nv.cid = newval.cid; 
00384   accessLock();
00385   theEdges[idx].update(ov, nv);
00386 #ifdef TDEBUG1
00387     CkPrintf("TMRC2D: [%d] Update edge %d replaced element %d with %d\n", cid,
00388          idx, ov.idx, nv.idx);
00389 #endif
00390   if (theEdges[idx].getBoundary() < b)
00391     theEdges[idx].setBoundary(b);
00392   if ((theEdges[idx].getBoundary() > 0) && (b > 0))
00393     CkPrintf("TMRC2D: [%d] WARNING! chunk::updateElement collapsed two different boundaries together on one edge.\n", cid);
00394   releaseLock();
00395 }
00396 
00397 void chunk::updateElementEdge(int idx, objRef oldval, objRef newval)
00398 {
00399   edgeRef ov, nv;
00400   ov.idx = oldval.idx;   ov.cid = oldval.cid; 
00401   nv.idx = newval.idx;   nv.cid = newval.cid; 
00402   accessLock();
00403   theElements[idx].update(ov, nv);
00404   releaseLock();
00405 }
00406 
00407 void chunk::updateReferences(int idx, objRef oldval, objRef newval)
00408 {
00409   CkPrintf("TMRC2D: [%d] WARNING! chunk::updateReferences called but not implemented!\n", cid);
00410   /*
00411   int i;
00412   nodeRef ov, nv;
00413   ov.idx = oldval.idx;   ov.cid = oldval.cid; 
00414   nv.idx = newval.idx;   nv.cid = newval.cid;
00415   for (i=0; i<numElements; i++)
00416     theElements[i].update(ov, nv);
00417   for (i=0; i<numEdges; i++)
00418     theEdges[i].updateSilent(ov, nv);
00419   */
00420 }
00421 
00422 doubleMsg *chunk::getArea(int n)
00423 {
00424   doubleMsg *dm = new doubleMsg;
00425   accessLock();
00426   dm->aDouble = theElements[n].getArea();
00427   releaseLock();
00428   return dm;
00429 }
00430 
00431 void chunk::resetEdge(int n)
00432 {
00433   accessLock();
00434   theEdges[n].reset();
00435   releaseLock();
00436 }
00437 
00438 refMsg *chunk::getNbr(int idx, objRef aRef)
00439 {
00440   refMsg *rm = new refMsg;
00441   elemRef er, ar;
00442   ar.cid = aRef.cid; ar.idx = aRef.idx;
00443   accessLock();
00444   er = theEdges[idx].getNot(ar);
00445   releaseLock();
00446   rm->aRef = er;
00447   return rm;
00448 }
00449 
00450 void chunk::setTargetArea(int idx, double aDouble)
00451 {
00452   accessLock();
00453   theElements[idx].setTargetArea(aDouble);
00454   releaseLock();
00455   if (!refineInProgress) {
00456     refineInProgress = 1;
00457     mesh[cid].refiningElements();
00458   }
00459 }
00460 
00461 void chunk::resetTargetArea(int idx, double aDouble)
00462 {
00463   accessLock();
00464   theElements[idx].resetTargetArea(aDouble);
00465   releaseLock();
00466 }
00467 
00468 void chunk::reportPos(int idx, double x, double y)
00469 {
00470   node z(x, y);
00471   accessLock();
00472   theNodes[idx].reportPos(z);
00473   releaseLock();
00474 }
00475 
00476 // the following methods are for run-time additions and modifications
00477 // to the chunk components
00478 void chunk::accessLock()
00479 {
00480   while (meshExpandFlag && (meshLock==0))
00481     CthYield();
00482   meshLock--;
00483 }
00484 
00485 void chunk::releaseLock()
00486 {
00487   meshLock++;
00488 }
00489 
00490 void chunk::adjustFlag()
00491 {
00492   meshExpandFlag = 1;
00493 }
00494 
00495 void chunk::adjustLock()
00496 {
00497   while (meshLock != 0)
00498     CthYield();
00499   meshLock = 1;
00500 }
00501 
00502 void chunk::adjustRelease()
00503 {
00504   meshLock = meshExpandFlag = 0;
00505 }
00506 
00507 void chunk::allocMesh(int nEl)
00508 {
00509   int i;
00510   sizeElements = nEl * 2;
00511   sizeNodes = sizeEdges = sizeElements * 3;
00512   elementSlots = nEl;
00513   firstFreeElement = nEl;
00514   theElements.resize(sizeElements);
00515   theNodes.resize(sizeNodes);
00516   theEdges.resize(sizeEdges);
00517   for (i=0; i<sizeElements; i++)
00518     theElements[i].set(); 
00519   for (i=0; i<sizeNodes; i++) {
00520     theNodes[i].present = 0; 
00521     theEdges[i].present = 0; 
00522   }
00523   conn = new int[3*numGhosts];
00524   gid = new int[2*numGhosts];
00525 }
00526 
00527 void chunk::adjustMesh()
00528 {
00529   int i;
00530   if (sizeElements <= elementSlots+100) {
00531     adjustFlag();
00532     adjustLock();
00533 #ifdef TDEBUG1
00534     CkPrintf("TMRC2D: [%d] Adjusting mesh size...\n", cid);
00535 #endif
00536     sizeElements += 100;
00537     sizeEdges += 300;
00538     sizeNodes += 300;
00539     theElements.resize(sizeElements);
00540     theEdges.resize(sizeEdges);
00541     theNodes.resize(sizeNodes);
00542     for (i=elementSlots; i<sizeElements; i++)
00543       theElements[i].present = 0;
00544     for (i=nodeSlots; i<sizeNodes; i++)
00545       theNodes[i].present = 0;
00546     for (i=edgeSlots; i<sizeEdges; i++)
00547       theEdges[i].present = 0;
00548 #ifdef TDEBUG1
00549     CkPrintf("TMRC2D: [%d] Done adjusting mesh size...\n", cid);
00550 #endif
00551     adjustRelease();
00552   }
00553 }
00554 
00555 intMsg *chunk::addNode(node n, int b1, int b2, int internal)
00556 {
00557   intMsg *im = new intMsg;
00558   im->anInt = firstFreeNode;
00559   theNodes[firstFreeNode] = n;
00560   theNodes[firstFreeNode].present = 1;
00561   if ((b1 == 0) || (b2 == 0) || internal) theNodes[firstFreeNode].boundary = 0;
00562   else if (b1 < b2) theNodes[firstFreeNode].boundary = b1; 
00563   else theNodes[firstFreeNode].boundary = b2;
00564   numNodes++;
00565   firstFreeNode++;
00566   if (firstFreeNode-1 == nodeSlots)  nodeSlots++;
00567   else  while (theNodes[firstFreeNode].isPresent()) firstFreeNode++;
00568   return im;
00569 }
00570 
00571 edgeRef chunk::addEdge(int n1, int n2, int b)
00572 {
00573 #ifdef TDEBUG1  
00574   CkPrintf("TMRC2D: [%d] Adding edge %d between nodes %d and %d\n", 
00575        cid, numEdges, n1, n2);
00576 #endif
00577   edgeRef eRef(cid, firstFreeEdge);
00578   theEdges[firstFreeEdge].set(firstFreeEdge, cid, this);
00579   theEdges[firstFreeEdge].reset();
00580   theEdges[firstFreeEdge].setNodes(n1, n2);
00581   theEdges[firstFreeEdge].setBoundary(b);
00582   numEdges++;
00583   firstFreeEdge++;
00584   if (firstFreeEdge-1 == edgeSlots)  edgeSlots++;
00585   else  while (theEdges[firstFreeEdge].isPresent()) firstFreeEdge++;
00586   return eRef;
00587 }
00588 
00589 edgeRef chunk::addEdge(int n1, int n2, const int *edgeBounds, const int *edgeConn, int nEdges)
00590 {
00591   int i=0;
00592 #ifdef TDEBUG1  
00593   CkPrintf("TMRC2D: [%d] Adding edge %d between nodes %d and %d\n", 
00594        cid, nEdges, n1, n2);
00595 #endif
00596   while (i<nEdges && !(n1==edgeConn[2*i] && n2==edgeConn[2*i+1]) || (n2==edgeConn[2*i] && n1==edgeConn[2*i+1]))
00597     i++;
00598   theEdges[i].set(i, cid, this);
00599   theEdges[i].reset();
00600   theEdges[i].setNodes(n1, n2);
00601   theEdges[i].setBoundary(edgeBounds[i]);
00602   edgeRef eRef(cid, i);
00603   numEdges++;
00604   while (theEdges[firstFreeEdge].isPresent()) firstFreeEdge++;
00605   return eRef;
00606 }
00607 
00608 
00609 elemRef chunk::addElement(int n1, int n2, int n3)
00610 {
00611   elemRef eRef(cid, firstFreeElement);
00612   theElements[firstFreeElement].set(cid, firstFreeElement, this);
00613   theElements[firstFreeElement].set(n1, n2, n3);
00614   theElements[firstFreeElement].calculateArea();
00615   numElements++;
00616   firstFreeElement++;
00617   if (firstFreeElement-1 == elementSlots)  elementSlots++;
00618   else  while (theElements[firstFreeElement].isPresent()) firstFreeElement++;
00619   if (!refineInProgress) {
00620     refineInProgress = 1;
00621     mesh[cid].refiningElements();
00622   }
00623   return eRef;
00624 }
00625 
00626 elemRef chunk::addElement(int n1, int n2, int n3,
00627               edgeRef er1, edgeRef er2, edgeRef er3)
00628 {
00629 #ifdef TDEBUG1
00630   CkPrintf("TMRC2D: [%d] New element %d added with nodes %d, %d and %d\n", 
00631        cid, firstFreeElement, n1, n2, n3);
00632 #endif
00633   elemRef eRef(cid, firstFreeElement);
00634   theElements[firstFreeElement].set(cid, firstFreeElement, this);
00635   theElements[firstFreeElement].set(n1, n2, n3, er1, er2, er3);
00636   theElements[firstFreeElement].calculateArea();
00637   numElements++;
00638   firstFreeElement++;
00639   if (firstFreeElement-1 == elementSlots)  elementSlots++;
00640   else  while (theElements[firstFreeElement].isPresent()) firstFreeElement++;
00641   if (!refineInProgress) {
00642     refineInProgress = 1;
00643     mesh[cid].refiningElements();
00644   }
00645   return eRef;
00646 }
00647 
00648 void chunk::removeNode(int n)
00649 {
00650 #ifdef TDEBUG1
00651   CkPrintf("TMRC2D: [%d] Removing node %d\n", cid, n);
00652 #endif
00653   if (theNodes[n].present) {
00654     theNodes[n].present = 0;
00655     theNodes[n].reset();
00656     if (n < firstFreeNode) firstFreeNode = n;
00657     else if ((n == firstFreeNode) && (firstFreeNode == nodeSlots)) {
00658       firstFreeNode--; nodeSlots--;
00659     }
00660     numNodes--;
00661   }
00662 }
00663 
00664 void chunk::removeEdge(int n)
00665 {
00666 #ifdef TDEBUG1
00667   CkPrintf("TMRC2D: [%d] Removing edge %d\n", cid, n);
00668 #endif
00669   if (theEdges[n].present) {
00670     theEdges[n].reset();
00671     theEdges[n].present = 0;
00672     theEdges[n].elements[0].reset();
00673     theEdges[n].elements[1].reset();
00674     theEdges[n].nodes[0] = theEdges[n].nodes[1] = -1;
00675     if (n < firstFreeEdge) firstFreeEdge = n;
00676     else if ((n == firstFreeEdge) && (firstFreeEdge == edgeSlots)) {
00677       firstFreeEdge--; edgeSlots--;
00678     }
00679     numEdges--;
00680   }
00681   else {
00682     CkPrintf("TMRC2D: [%d] WARNING: chunk::removeEdge(%d): edge not present\n", cid, n);
00683   } 
00684 }
00685 
00686 void chunk::removeElement(int n)
00687 {
00688 #ifdef TDEBUG1
00689   CkPrintf("TMRC2D: [%d] Removing element %d\n", cid, n);
00690 #endif
00691   if (theElements[n].present) {
00692     theElements[n].present = 0;
00693     theElements[n].clear();
00694     if (n < firstFreeElement) firstFreeElement = n;
00695     else if ((n == firstFreeElement) && (firstFreeElement == elementSlots)) {
00696       firstFreeElement--; elementSlots--;
00697     }
00698     numElements--;
00699   }
00700   else {
00701     CkPrintf("TMRC2D: [%d] WARNING: chunk::removeElement(%d): element not present\n", cid, n);
00702   } 
00703 }
00704 
00705 // these two functions produce debugging output by printing somewhat
00706 // sychronized versions of the entire mesh to files readable by tkplotter
00707 void chunk::print()
00708 {
00709   accessLock();
00710   debug_print(debug_counter);
00711   debug_counter++;
00712   releaseLock();
00713 }
00714 
00715 void chunk::debug_print(int c)
00716 {
00717   FILE *fp;
00718   char filename[30];
00719   int i, j;
00720   node n;
00721 
00722   memset(filename, 0, 30);
00723   sprintf(filename, "mesh_debug_%d.%d", cid, c);
00724   fp = fopen(filename, "w");
00725 
00726   fprintf(fp, "%d %d\n", cid, numElements);
00727   for (i=0; i<elementSlots; i++) {
00728     if (theElements[i].isPresent()) {
00729       for (j=0; j<3; j++) {
00730     n = theNodes[theElements[i].getNode(j)];
00731     fprintf(fp, "%f %f   ", n.X(), n.Y());
00732       }
00733       fprintf(fp, "%d %f\n", i, theElements[i].getTargetArea());
00734       for (j=0; j<3; j++) {
00735     fprintf(fp, "0  ");
00736       }
00737       fprintf(fp, "\n");
00738     }
00739   }
00740   fclose(fp);
00741 }
00742 
00743 void chunk::out_print()
00744 {
00745   FILE *fp;
00746   char filename[30];
00747   int i;
00748 
00749   memset(filename, 0, 30);
00750   sprintf(filename, "mesh.out");
00751   fp = fopen(filename, "a");
00752 
00753   if (cid == 0)
00754     fprintf(fp, "%d\n", numChunks);
00755   fprintf(fp, " %d %d %d %d\n", cid, numNodes, numEdges, numElements);
00756   for (i=0; i<numNodes; i++)
00757     fprintf(fp, "    %f %f\n", theNodes[i].X(), theNodes[i].Y());
00758   for (i=0; i<numEdges; i++) {
00759     fprintf(fp, " %d %d ", theEdges[i].elements[0].idx, theEdges[i].elements[0].cid);
00760     fprintf(fp, " %d %d\n", theEdges[i].elements[1].idx, theEdges[i].elements[1].cid);
00761   }
00762   for (i=0; i<numElements; i++) {
00763     if (theElements[i].isPresent()) {
00764       fprintf(fp, " %d ", theElements[i].nodes[0]);
00765       fprintf(fp, " %d ", theElements[i].nodes[1]);
00766       fprintf(fp, " %d ", theElements[i].nodes[2]);
00767       fprintf(fp, "   ");
00768       fprintf(fp, " %d %d ", theElements[i].edges[0].idx, theElements[i].edges[0].cid);
00769       fprintf(fp, " %d %d ", theElements[i].edges[1].idx, theElements[i].edges[1].cid);
00770       fprintf(fp, " %d %d\n", theElements[i].edges[2].idx, theElements[i].edges[2].cid);
00771     }
00772   }
00773   fprintf(fp, "\n");
00774   fclose(fp);
00775 }
00776 
00777 void chunk::updateNodeCoords(int nNode, double *coord, int nEl)
00778 {
00779   int i;
00780 #ifdef TDEBUG2
00781   CkPrintf("TMRC2D: [%d] updateNodeCoords...\n", cid);
00782 #endif
00783   if (first) {
00784     CkWaitQD();
00785     first = 0;
00786   }
00787 #ifdef TDEBUG3
00788   CkPrintf("TMRC2D: [%d] In updateNodeCoords after CkWaitQD: edges sent=%d edge recvd=%d\n", cid, edgesSent, edgesRecvd);
00789 #endif
00790 #ifdef TDEBUG1
00791   sanityCheck(); // quietly make sure mesh is in shape
00792   validCheck();
00793 #endif
00794   // do some error checking
00795   if (nEl != numElements) {
00796     CkPrintf("TMRC2D: [%d] ERROR: nEl=%d passed in updateNodeCoords does not match TMRC2D numElements=%d!\n", cid, nEl, numElements);
00797     CkAssert(nEl == numElements);
00798   }
00799   if (nNode != numNodes) {
00800     CkPrintf("TMRC2D: [%d] ERROR: nNode=%d passed in updateNodeCoords does not match TMRC2D numNodes=%d!\n", cid, nNode, numNodes);
00801     CkAssert(nNode == numNodes);
00802   } 
00803   // update node coordinates from coord
00804   for (i=0; i<nodeSlots; i++)
00805     if (theNodes[i].isPresent()) {
00806       theNodes[i].set(coord[2*i], coord[2*i + 1]);
00807 #ifdef TDEBUG3      
00808       if (theNodes[i].boundary) {
00809     CkPrintf("TMRC2D: [%d] Node %d on boundary!\n", cid, i);
00810       } 
00811 #endif
00812     }
00813   // recalculate and cache new areas for each element
00814   for (i=0; i<elementSlots; i++) 
00815     if (theElements[i].isPresent())  theElements[i].calculateArea();
00816 #ifdef TDEBUG1
00817   sanityCheck(); // quietly make sure mesh is in shape
00818   dump();
00819 #endif
00820 #ifdef TDEBUG2
00821   CkPrintf("TMRC2D: [%d] updateNodeCoords DONE.\n", cid);
00822 #endif
00823 }
00824 
00825 void chunk::multipleRefine(double *desiredArea, refineClient *client)
00826 {
00827   int i;
00828   double area;
00829 #ifdef TDEBUG2
00830   CkPrintf("TMRC2D: [%d] multipleRefine....\n", cid);
00831 #endif
00832 #ifdef TDEBUG1
00833   sanityCheck(); // quietly make sure mesh is in shape
00834 #endif
00835   CmiMemoryCheck();
00836   theClient = client; // initialize refine client associated with this chunk
00837   
00838   if (refineStack) delete [] refineStack;
00839   refineStack=new elemHeap[numElements];
00840   if (refineElements) delete [] refineElements;
00841   refineElements = new elemHeap[numElements+1];
00842   for (i=0; i<elementSlots; i++) { // set desired areas for elements
00843     if (theElements[i].isPresent()) {
00844       area = theElements[i].getArea();
00845       if (area == 0.0) {
00846     CkPrintf("Element[%d] has area %1.10e!\n", i, area);
00847       }
00848       //precThrshld = area * 1e-8;
00849       if (desiredArea[i] < REFINE_PRECISION*area) {
00850     theElements[i].resetTargetArea(desiredArea[i]);
00851     double qFactor=theElements[i].getAreaQuality();
00852     Insert(i, qFactor, 0);
00853     CkPrintf("Element[%d] has area %1.10e target %1.10e qFactor %1.10e\n", i, area, desiredArea[i], qFactor);
00854 #ifdef TDEBUG2
00855     CkPrintf("TMRC2D: [%d] Setting target on element %d to %1.10e with largeEdge %1.10e\n", cid, i, desiredArea[i], qFactor);
00856 #endif
00857       }
00858     }
00859   }
00860   // start refinement
00861   if (!refineInProgress) {
00862    refineInProgress = 1;
00863    mesh[cid].refiningElements();
00864   }
00865 #ifdef TDEBUG2
00866   CkPrintf("TMRC2D: [%d] multipleRefine DONE.\n", cid);
00867 #endif
00868   CkWaitQD();
00869   CmiMemoryCheck();
00870 }
00871 
00872 void chunk::multipleCoarsen(double *desiredArea, refineClient *client)
00873 {
00874   int i;
00875   double precThrshld, area, qFactor;
00876   CmiMemoryCheck();
00877 #ifdef TDEBUG2
00878   CkPrintf("TMRC2D: [%d] multipleCoarsen....\n", cid);
00879 #endif
00880 #ifdef TDEBUG1
00881   sanityCheck(); // quietly make sure mesh is in shape
00882 #endif
00883   theClient = client; // initialize refine client associated with this chunk
00884   if (coarsenElements) delete [] coarsenElements;
00885   coarsenElements = new elemHeap[numElements+1];
00886   coarsenElements[0].len=-2.0;
00887   coarsenElements[0].elID=-1;
00888   for (i=0; i<elementSlots; i++) { // set desired areas for elements
00889     if (theElements[i].isPresent()) {
00890       area = theElements[i].getArea();
00891       if (area == 0.0) {
00892     CkPrintf("Element[%d] has area %1.10e!\n", i, area);
00893       }
00894       double shortEdgeLen, angle;
00895       shortEdgeLen = theElements[i].getShortestEdge(&angle);
00896       qFactor = theElements[i].getAreaQuality();
00897       theElements[i].nonCoarsenCount = 0;
00898       //precThrshld = area * 1e-8;
00899       theElements[i].resetTargetArea(desiredArea[i]);
00900       if ((desiredArea[i]*COARSEN_PRECISION > area) || (area == 0.0) || 
00901       (qFactor<QUALITY_MIN)) {
00902     Insert(i, shortEdgeLen*qFactor, 1);
00903     CkPrintf("Element[%d] has area %1.10e target %1.10e qFactor %1.10e\n", i, area, desiredArea[i], qFactor);
00904 #ifdef TDEBUG2
00905     CkPrintf("TMRC2D: [%d] Setting target on element %d to %1.10e with shortEdge %1.10e\n", cid, i, desiredArea[i], qFactor);
00906 #endif
00907       }
00908     }
00909   }
00910   
00911   // start coarsening
00912   if (!coarsenInProgress) {
00913    coarsenInProgress = 1;
00914    mesh[cid].coarseningElements();
00915   }
00916 #ifdef TDEBUG2
00917   CkPrintf("TMRC2D: [%d] multipleCoarsen DONE.\n", cid);
00918 #endif
00919   CkWaitQD();
00920 /*
00921   for (i=0; i<elementSlots; i++) { // check desired areas for elements
00922     if (theElements[i].isPresent()) {
00923       area = theElements[i].getArea();
00924       if (desiredArea[i] > area) {
00925     CkPrintf("TMRC2D: [%d] WARNING: element %d area is %1.10e but target was %1.10e\n", cid, i, area, desiredArea[i]);
00926       }
00927     }
00928   }
00929 */
00930   CmiMemoryCheck();
00931 }
00932 
00933 void chunk::newMesh(int meshID_,int nEl, int nGhost, const int *conn_, 
00934             const int *gid_, int nnodes, const int *boundaries, int nEdges, 
00935             const int *edgeConn, const int *edgeBounds, int idxOffset)
00936 {
00937   int i, j;
00938 #ifdef TDEBUG2
00939   CkPrintf("TMRC2D: [%d] In newMesh...\n", cid);
00940 #endif
00941   meshID = meshID_;
00942 
00943   meshPtr = FEM_Mesh_lookup(meshID, "chunk::newMesh");
00944   numElements=nEl;
00945   numGhosts = nGhost;
00946   allocMesh(nEl);
00947   int *conn = new int[3*numGhosts];
00948   int *gid = new int[2*numGhosts];
00949 
00950   // add elements to chunk
00951   for (i=0; i<numElements; i++) {
00952     int nodes[3];
00953     edgeRef edges[3];
00954      for (j=0; j<3; j++) {
00955       int c=conn_[i*3+j]-idxOffset;
00956       conn[i*3 + j] = c;
00957       nodes[j] = c;
00958       edges[j].reset();
00959      }
00960     theElements[i].set(cid, i, this);
00961     theElements[i].set(nodes, edges);
00962     gid[i*2] = cid;
00963     gid[i*2 + 1] = i;
00964   }
00965 
00966   // add ghost elements to chunk
00967   for (i=nEl; i<nGhost; i++) {
00968     for (j=0; j<3; j++)
00969       conn[i*3+j] = conn_[i*3+j]-idxOffset;
00970     gid[i*2+0] = gid_[i*2]-idxOffset;
00971     gid[i*2+1] = gid_[i*2 + 1]-idxOffset;
00972   }
00973 
00974   MPI_Barrier(MPI_COMM_WORLD);
00975   // derive edges from elements on this chunk
00976   if (boundaries) {
00977     numNodes = nnodes;
00978 #ifdef TDEBUG1
00979     CkPrintf("TMRC2D: [%d] Received non-NULL boundaries... adding...\n", cid);
00980 #endif
00981     for (i=0; i<numNodes; i++) {
00982       theNodes[i].boundary = boundaries[i];
00983 #ifdef TDEBUG3
00984       if (theNodes[i].boundary) 
00985     CkPrintf("TMRC2D: [%d] Node %d has boundary %d.\n", cid, i, theNodes[i].boundary);
00986 #endif
00987     }
00988     deriveEdges(conn, gid, edgeBounds, edgeConn, nEdges);
00989     CkAssert(nnodes == numNodes);
00990   }
00991   else {
00992     deriveEdges(conn, gid, edgeBounds, edgeConn, nEdges);
00993     CkAssert(nnodes == numNodes);
00994     deriveBoundaries(conn, gid);
00995   }
00996   delete[] conn;
00997   delete[] gid;
00998 #ifdef TDEBUG2
00999   CkPrintf("TMRC2D: [%d] newMesh DONE; chunk created with %d elements.\n", 
01000        cid, numElements);
01001 #endif
01002 }
01003 
01004 void chunk::deriveEdges(int *conn, int *gid, const int *edgeBounds, const int *edgeConn, int nEdges)
01005 {
01006   // need to add edges to the chunk, and update all edgeRefs on all elements
01007   // also need to add nodes to the chunk
01008   int i, j, n1localIdx, n2localIdx;
01009   edgeRef newEdge;
01010 
01011   deriveNodes(); // now numNodes and theNodes have values
01012 #ifdef TDEBUG2
01013   CkPrintf("TMRC2D: [%d] Deriving edges...\n", cid);
01014 #endif
01015   for (i=0; i<elementSlots; i++) {
01016 #ifdef TDEBUG3
01017     CkPrintf("TMRC2D: [%d] Deriving edges for element %d...\n", cid, i);
01018 #endif
01019     elemRef myRef(cid,i);
01020     for (j=0; j<3; j++) {
01021       n1localIdx = j;
01022       n2localIdx = (j+1) % 3;
01023 #ifdef TDEBUG3
01024     CkPrintf("TMRC2D: [%d] Deriving edges for element %d between %d,%d (real nodes %d,%d)...\n", cid, i, n1localIdx, n2localIdx, theElements[i].nodes[n1localIdx], theElements[i].nodes[n2localIdx]);
01025 #endif
01026       // look for edge
01027       if (theElements[i].edges[j] == nullRef) { // the edge doesn't exist yet
01028     // get nbr ref
01029 #ifdef TDEBUG3
01030     CkPrintf("TMRC2D: [%d] Edge between nodes %d,%d doesn't exist yet...\n", cid, theElements[i].nodes[n1localIdx], theElements[i].nodes[n2localIdx]);
01031 #endif
01032     elemRef nbrRef;
01033     int edgeIdx = getNbrRefOnEdge(theElements[i].nodes[n1localIdx], 
01034                       theElements[i].nodes[n2localIdx], 
01035                       conn, numGhosts, gid, i, &nbrRef);
01036     FEM_Node *mNodes = &(meshPtr->node);
01037     int nIdx1 = theElements[i].nodes[n1localIdx];
01038     int nIdx2 = theElements[i].nodes[n2localIdx];
01039     if ((theNodes[nIdx1].boundary < theNodes[nIdx2].boundary) && 
01040         (theNodes[nIdx1].boundary != 0))    {
01041       theNodes[nIdx2].fixed = 1;
01042 #ifdef TDEBUG2
01043       CkPrintf("TMRC2D: [%d] Corner detected, fixing node %d\n", 
01044            cid, nIdx2);
01045 #endif
01046       FEM_Comm_Rec *nRec = (FEM_Comm_Rec*)(mNodes->shared.getRec(nIdx2));
01047       if (nRec) {
01048         int count = nRec->getShared();
01049         for (i=0; i<count; i++) {
01050           mesh[nRec->getChk(i)].fixNode(nRec->getIdx(i), cid);
01051         }
01052       }
01053     }
01054     if ((theNodes[nIdx2].boundary < theNodes[nIdx1].boundary) && 
01055         (theNodes[nIdx2].boundary != 0)) {
01056       theNodes[nIdx1].fixed = 1;
01057 #ifdef TDEBUG2
01058       CkPrintf("TMRC2D: [%d] Corner detected, fixing node %d\n", 
01059            cid, nIdx1);
01060 #endif
01061       FEM_Comm_Rec *nRec = (FEM_Comm_Rec*)(mNodes->shared.getRec(nIdx1));
01062       if (nRec) {
01063         int count = nRec->getShared();
01064         for (i=0; i<count; i++) {
01065           mesh[nRec->getChk(i)].fixNode(nRec->getIdx(i), cid);
01066         }
01067       }
01068     }
01069     if (edgeLocal(myRef, nbrRef)) { // make edge here
01070       if (!edgeBounds)
01071         newEdge = addEdge(nIdx1, nIdx2, 0);
01072       else
01073         newEdge = addEdge(nIdx1, nIdx2, edgeBounds, edgeConn, nEdges);
01074 #ifdef TDEBUG3
01075       CkPrintf("TMRC2D: [%d] New edge (%d,%d) added between nodes %d and %d and elements %d and %d on chunk %d\n", cid, newEdge.cid, newEdge.idx, nIdx1, nIdx2, i, nbrRef.idx, nbrRef.cid);
01076 #endif
01077       // point edge to the two neighboring elements
01078       theEdges[newEdge.idx].update(nullRef, myRef);
01079       theEdges[newEdge.idx].update(nullRef, nbrRef);
01080       // point elem i's edge j at the edge
01081       theElements[i].set(j, newEdge);
01082       // point nbrRef at the edge
01083       if (nbrRef.cid==cid) // Local neighbor
01084         theElements[nbrRef.idx].set(edgeIdx, newEdge);
01085       else if (nbrRef.cid != -1) { // Remote neighbor
01086         mesh[nbrRef.cid].addRemoteEdge(nbrRef.idx, edgeIdx, newEdge);
01087         edgesSent++;
01088       }
01089     }
01090         else { // else edge will be made on a different chunk
01091 #ifdef TDEBUG3
01092       CkPrintf("TMRC2D: [%d] Edge between nodes %d,%d to be created elsewhere...\n", cid, theElements[i].nodes[n1localIdx], theElements[i].nodes[n2localIdx]);
01093 #endif
01094       // mark the edge as non-boundary; will be filled later
01095       theElements[i].edges[j].idx = theElements[i].edges[j].cid = -2;
01096     }
01097       }
01098 #ifdef TDEBUG3
01099       else { 
01100     CkPrintf("TMRC2D: [%d] Edge between nodes %d,%d exists at %d\n", cid, theElements[i].nodes[n1localIdx], theElements[i].nodes[n2localIdx], theElements[i].edges[j].idx);
01101       }
01102 #endif
01103     }
01104   }
01105   nodeSlots = numNodes;
01106   firstFreeNode = numNodes;
01107   edgeSlots = numEdges;
01108   firstFreeEdge = numEdges;
01109 #ifdef TDEBUG2
01110   CkPrintf("TMRC2D: [%d] Done deriving edges...\n", cid);
01111 #endif
01112 }
01113 
01114 void chunk::deriveNodes()
01115 {
01116   int i, j;
01117   int aNode;
01118 
01119 #ifdef TDEBUG2
01120   CkPrintf("TMRC2D: [%d] Deriving nodes...\n", cid);
01121 #endif
01122   numNodes = 0;
01123   for (i=0; i<elementSlots; i++) {
01124     for (j=0; j<3; j++) {
01125       aNode = theElements[i].nodes[j];
01126       CkAssert(aNode > -1);
01127       if ((aNode + 1) > numNodes)
01128     numNodes = aNode + 1;
01129       theNodes[aNode].present = 1;
01130     }
01131   }
01132 #ifdef TDEBUG3
01133   CkPrintf("TMRC2D: [%d] NumNodes = %d; max node idx = %d\n", 
01134        cid, numNodes, numNodes-1);
01135 #endif
01136 #ifdef TDEBUG2
01137   CkPrintf("TMRC2D: [%d] Done deriving nodes.\n", cid);
01138 #endif
01139 }
01140 
01141 int chunk::edgeLocal(elemRef e1, elemRef e2)
01142 {
01143   return ((e1.cid==-1 || e2.cid==-1) || (e1.cid > e2.cid) ||
01144       ((e1.cid == e2.cid) && (e1.idx < e2.idx)));
01145 }
01146 
01147 int chunk::getNbrRefOnEdge(int n1, int n2, int *conn, int nGhost, int *gid, 
01148                int idx, elemRef *er)
01149 {
01150   int i, e;
01151   er->set(-1, -1);
01152   for (i=0; i<nGhost; i++) {
01153     if (i != idx) {
01154       if ((e = hasEdge(n1, n2, conn, i)) != -1) {
01155     er->set(gid[i*2], gid[i*2+1]);
01156     return e;
01157       }
01158     }
01159   }
01160   return -1;
01161 }
01162 
01163 int chunk::hasEdge(int n1, int n2, int *conn, int idx) 
01164 {
01165   int i, j;
01166   for (i=0; i<3; i++) {
01167     j = (i+1) % 3;
01168     if (((conn[idx*3+i] == n1) && (conn[idx*3+j] == n2)) ||
01169     ((conn[idx*3+j] == n1) && (conn[idx*3+i] == n2)))
01170       return i;
01171   }
01172   return -1;
01173 }
01174 
01175 void chunk::deriveBoundaries(int *conn, int *gid)
01176 {
01177   int edgeIdx;
01178   elemRef nbrRef;
01179   CkPrintf("TMRC2D: [%d] WARNING! Null list of boundary flags passed to newMesh...\n ...I hope you didn't want coarsening to work!\n", cid);
01180   CkPrintf("TMRC2D: [%d] WARNING! Attempting to derive single boundary information -- corners will be ignored!\n", cid);
01181   for (int i=0; i<numElements; i++) {
01182     for (int j=0; j<3; j++) {
01183       edgeIdx = getNbrRefOnEdge(theElements[i].nodes[j], 
01184                 theElements[i].nodes[(j+1)%3], 
01185                 conn, numGhosts, gid, i, &nbrRef); 
01186       if ((nbrRef.idx == -1) && (nbrRef.cid == -1)) {
01187 #ifdef TDEBUG2
01188     CkPrintf("TMRC2D: [%d] Marking node %d as a boundary.\n", 
01189          cid, theElements[i].nodes[j]);
01190     CkPrintf("TMRC2D: [%d] Marking node %d as a boundary.\n", 
01191          cid, theElements[i].nodes[(j+1)%3]);
01192 #endif
01193     theNodes[theElements[i].nodes[j]].boundary = 1;
01194     theNodes[theElements[i].nodes[(j+1)%3]].boundary = 1;
01195       }
01196     }
01197   }
01198 }
01199 
01200 void chunk::tweakMesh()
01201 {
01202   CkPrintf("TMRC2D: [%d] WARNING! chunk::tweakMesh called but not implemented!\n", cid);
01203   /*
01204   for (int i=0; i<numElements; i++) 
01205     theElements[i].tweakNodes();
01206   */
01207 }
01208 
01209 void chunk::improveChunk()
01210 {
01211   CkPrintf("TMRC2D: [%d] WARNING! chunk::improveChunk called but not implemented!\n", cid);
01212   /*
01213   for (int i=0; i<numNodes; i++) 
01214     if (!theNodes[i].border) {
01215       theNodes[i].improvePos();
01216     }
01217     else CkPrintf("Not adjusting node %d on chunk %d\n", i, thisIndex);
01218   */
01219 }
01220 
01221 void chunk::improve()
01222 {
01223   CkPrintf("TMRC2D: [%d] WARNING! chunk::improve called but not implemented!\n", cid);
01224   /*
01225   for (int i=0; i<20; i++) {
01226     mesh.tweakMesh();
01227     mesh.improveChunk();
01228   }
01229   */
01230 }
01231 
01232 void chunk::sanityCheck(void)
01233 {
01234   CkPrintf("TMRC2D: [%d] running sanity check...\n", cid);
01235   int i;
01236   if (numElements<0 || (int)theElements.size()<numElements)
01237         CkAbort("sc-> TMRC2D: numElements or vector size insane!");
01238   if (numEdges<0 || (int)theEdges.size()<numEdges)
01239         CkAbort("sc-> TMRC2D: numEdges or vector size insane!");
01240   if (numNodes<0 || (int)theNodes.size()<numNodes)
01241         CkAbort("sc-> TMRC2D: numNodes or vector size insane!");
01242   i=0; 
01243   while (theElements[i].isPresent()) i++;
01244   CkAssert(firstFreeElement == i);
01245   i=0; 
01246   while (theEdges[i].isPresent()) i++;
01247   CkAssert(firstFreeEdge == i);
01248   i=0; 
01249   while (theNodes[i].isPresent()) i++;
01250   CkAssert(firstFreeNode == i);
01251   for (i=0; i<elementSlots; i++) 
01252     if (theElements[i].isPresent())
01253       theElements[i].sanityCheck(this, elemRef(cid,i), nodeSlots);
01254   for (i=0; i<edgeSlots; i++)
01255     if (theEdges[i].isPresent())
01256       theEdges[i].sanityCheck(this, edgeRef(cid,i));
01257   for (i=0; i<nodeSlots; i++)
01258     if (theNodes[i].isPresent())
01259       theNodes[i].sanityCheck(cid, i);
01260   CkPrintf("TMRC2D: [%d] sanity check PASSED.\n", cid);
01261 }
01262 
01263 void chunk::dump()
01264 {
01265   int i;
01266   CkPrintf("TMRC2D: [%d] has %d elements, %d nodes and %d edges:\n", 
01267        cid, numElements, numNodes, numEdges);
01268   for (i=0; i<elementSlots; i++) 
01269     if (theElements[i].isPresent()){
01270       CkPrintf("TMRC2D: [%d] Element %d nodes %d %d %d edges [%d,%d] [%d,%d] [%d,%d]\n", cid, i, theElements[i].nodes[0], theElements[i].nodes[1], theElements[i].nodes[2], theElements[i].edges[0].cid, theElements[i].edges[0].idx, theElements[i].edges[1].cid, theElements[i].edges[1].idx, theElements[i].edges[2].cid, theElements[i].edges[2].idx);
01271         }   
01272   for (i=0; i<nodeSlots; i++)
01273     if (theNodes[i].isPresent()) {
01274       CkPrintf("TMRC2D: [%d] Node %d has coordinates (%f,%f)\n", 
01275            cid, i, theNodes[i].X(), theNodes[i].Y());
01276     }        
01277   for (i=0; i<edgeSlots; i++)
01278     if (theEdges[i].isPresent()) {
01279       CkPrintf("TMRC2D: [%d] Edge %d has nodes %d %d elements [%d,%d] [%d,%d]\n", cid, 
01280            i, theEdges[i].nodes[0], theEdges[i].nodes[1], theEdges[i].elements[0].cid, theEdges[i].elements[0].idx, 
01281            theEdges[i].elements[1].cid, theEdges[i].elements[1].idx);
01282     }        
01283 }
01284 
01285 intMsg *chunk::lockChunk(int lhc, int lhi, double prio)
01286 {
01287   intMsg *im = new intMsg;
01288   im->anInt = lockLocalChunk(lhc, lhi, prio);
01289   return im;
01290 }
01291 
01292 int chunk::lockLocalChunk(int lhc, int lhi, double prio)
01293 { // priority is given to the lower prio, lower lh
01294   // each chunk can have at most one lock reservation
01295   // if a chunk holds a lock, it cannot have a reservation
01296   if (!lock) { // this chunk is UNLOCKED
01297     // remove previous lh reservation; insert this reservation
01298     // then check first reservation
01299     removeLock(lhc, lhi);
01300     insertLock(lhc, lhi, prio);
01301     if ((lockList->prio == prio) && (lockList->holderIdx == lhi)
01302     && (lockList->holderCid == lhc)) {
01303 #ifdef TDEBUG3
01304       CkPrintf("TMRC2D: [%d] LOCK chunk %d by %d on %d prio %.10f LOCKED(was %d[%don%d:%.10f])\n", CkMyPe(), cid, lhi, lhc, prio, lock, lockHolderIdx, lockHolderCid, lockPrio);
01305 #endif
01306       lock = 1; lockHolderCid = lhc; lockHolderIdx = lhi; lockPrio = prio; 
01307       //lockCount = 1;
01308       lockList = lockList->next;
01309       return 1;
01310     }
01311     removeLock(lhc, lhi);
01312 #ifdef TDEBUG3
01313     CkPrintf("TMRC2D: [%d] LOCK chunk %d by %d on %d prio %.10f REFUSED (was %d[%don%d:%.10f])\n", CkMyPe(), cid, lhi, lhc, prio, lock, lockHolderIdx, lockHolderCid, lockPrio);
01314 #endif
01315     return 0;
01316   }
01317   else { // this chunk is LOCKED
01318     if ((prio == lockPrio) && (lhi == lockHolderIdx) && (lhc == lockHolderCid)) {
01319       //lockCount++;
01320 #ifdef TDEBUG3
01321       CkPrintf("TMRC2D: [%d] LOCK chunk %d by %d on %d prio %.10f HELD (was %d[%don%d:%.10f])\n", CkMyPe(), cid, lhi, lhc, prio, lock, lockHolderIdx, lockHolderCid, lockPrio);
01322 #endif
01323       return 1;
01324     }
01325     if ((lhi == lockHolderIdx) && (lhc == lockHolderCid)) {
01326       CkPrintf("ERROR: lockholder %d on %d trying to relock with different prio! prio=%f newprio=%f\n", lhi, lhc, lockPrio, prio);
01327       CmiAssert(!((lhi == lockHolderIdx) && (lhc == lockHolderCid)));
01328     }
01329     /*
01330     removeLock(lhc, lhi);
01331     insertLock(lhc, lhi, prio);
01332     if ((lockList->prio == prio) && (lockList->holderIdx == lhi) && (lockList->holderCid == lhc)) {
01333       if ((prio < lockPrio) && (lhc != lockHolderCid)) { // lh might be next
01334     CkPrintf("[%d]LOCK chunk %d by %d on %d prio %.10f SPIN (was %d[%don%d:%.10f])\n", CkMyPe(), cid, lhi, lhc, prio, lock, lockHolderIdx, lockHolderCid, lockPrio);
01335     return 0; 
01336       }
01337     }
01338     removeLock(lhc, lhi);
01339     */
01340 #ifdef TDEBUG3
01341     CkPrintf("TMRC2D: [%d] LOCK chunk %d by %d on %d prio %.10f REFUSED (was %d[%don%d:%.10f])\n", CkMyPe(), cid, lhi, lhc, prio, lock, lockHolderIdx, lockHolderCid, lockPrio);
01342 #endif
01343     return 0;
01344   }
01345 }
01346 
01347 void chunk::insertLock(int lhc, int lhi, double prio)
01348 {
01349   prioLockStruct *newLock, *tmp;
01350   
01351   newLock = (prioLockStruct *)malloc(sizeof(prioLockStruct));
01352   newLock->prio = prio;
01353   newLock->holderIdx = lhi;
01354   newLock->holderCid = lhc;
01355   newLock->next = NULL;
01356   if (!lockList) lockList = newLock;
01357   else {
01358     if ((prio < lockList->prio) || 
01359     ((prio == lockList->prio) && (lhc < lockList->holderCid)) ||
01360     ((prio == lockList->prio) && (lhc == lockList->holderCid) && 
01361      (lhi < lockList->holderIdx))) {
01362       // insert before first element in lockList
01363       newLock->next = lockList;
01364       lockList = newLock;
01365     }
01366     else { 
01367       tmp = lockList;
01368       while (tmp->next) {
01369     if ((prio > tmp->next->prio) || 
01370         ((prio == tmp->next->prio) && (lhc > tmp->next->holderCid)) ||
01371         ((prio == tmp->next->prio) && (lhc == tmp->next->holderCid) &&
01372          (lhi > tmp->next->holderIdx))) 
01373       tmp = tmp->next;
01374     else break;
01375       }
01376       newLock->next = tmp->next;
01377       tmp->next = newLock;
01378     }
01379   }
01380 }
01381 
01382 void chunk::removeLock(int lhc, int lhi) 
01383 {
01384   prioLockStruct *tmp = lockList, *del;
01385   if (!lockList) return;
01386   if ((tmp->holderCid == lhc) && (tmp->holderIdx == lhi)) {
01387     lockList = lockList->next;
01388     free(tmp);
01389   }
01390   else {
01391     while (tmp->next && 
01392        ((tmp->next->holderCid != lhc) || (tmp->next->holderIdx != lhi)))
01393       tmp = tmp->next;
01394     if (tmp->next) {
01395       del = tmp->next;
01396       tmp->next = del->next;
01397       free(del);
01398     }
01399   }
01400 }
01401 
01402 void chunk::unlockChunk(int lhc, int lhi)
01403 {
01404   unlockLocalChunk(lhc, lhi);
01405 }
01406 
01407 void chunk::unlockLocalChunk(int lhc, int lhi)
01408 {
01409   if ((lockHolderIdx == lhi) && (lockHolderCid == lhc)) {
01410     //    lockCount--;
01411     //    if (lockCount == 0) {
01412       lock = 0; 
01413       lockHolderIdx = -1;
01414       lockHolderCid = -1;
01415       lockPrio = -1.0;
01416 #ifdef TDEBUG3
01417       CkPrintf("TMRC2D: [%d] UNLOCK chunk %d by holder %d on %d\n", 
01418            CkMyPe(), cid, lhi, lhc);
01419 #endif
01420     //    }
01421   }
01422   /*  else { 
01423     CkPrintf("TMRC2D: [%d] ERROR: UNLOCK chunk %d by holder %d on %d -- THIS WAS NOT THE HOLDER!!! %d on %d was!\n", CkMyPe(), cid, lhi, lhc, lockHolderIdx, lockHolderCid);
01424     CkAbort(0);
01425     }*/
01426 }
01427 
01428 void chunk::fixNode(int nIdx, int chkid)
01429 {
01430   FEM_Node *nodesPtr = &(meshPtr->node);
01431   const FEM_Comm_List *sharedList = &(nodesPtr->shared.getList(chkid));
01432   nIdx = (*sharedList)[nIdx];
01433 
01434   theNodes[nIdx].fixed = 1;
01435 #ifdef TDEBUG2
01436   CkPrintf("TMRC2D: [%d] Corner detected, fixing node %d\n", cid, nIdx);
01437 #endif
01438 }
01439 
01440 int chunk::joinCommLists(int nIdx, int shd, int *chk, int *idx, int *rChk,
01441              int *rIdx)
01442 {
01443   FEM_Node *theNodes = &(meshPtr->node);
01444   FEM_Comm_Rec *nodeRec=(FEM_Comm_Rec *)(theNodes->shared.getRec(nIdx));
01445   int nShared, count, chunk, index, found, i, j;
01446   if (nodeRec) nShared = nodeRec->getShared();
01447   else { 
01448     rChk = chk;
01449     rIdx = idx;
01450     return shd;
01451   }
01452   rChk = (int *)malloc((shd+nShared)*sizeof(int));
01453   rIdx = (int *)malloc((shd+nShared)*sizeof(int));
01454   for (i=0; i<shd; i++) {
01455     rChk[i] = chk[i];
01456     rIdx[i] = idx[i];
01457   }
01458   count = shd;
01459   for (i=0; i<nShared; i++) {
01460     chunk = nodeRec->getChk(i);
01461     for (j=0; j<shd; j++) {
01462       if (rChk[j] == chunk) { // already present
01463     found = 1;
01464     break;
01465       }
01466     }
01467     if (!found) {
01468       index = nodeRec->getIdx(i-shd);
01469       rChk[count] = chunk;
01470       rIdx[count] = index;
01471       count++;
01472     }
01473   }
01474   free(rChk);
01475   free(rIdx);
01476   return count;
01477 }
01478 
01479 void chunk::Insert(int eIdx, double len, int cFlag)
01480 {
01481   int i;
01482   if (cFlag) {
01483     i = ++coarsenHeapSize; 
01484     while ((coarsenElements[i/2].len>=len) && (i != 1))
01485     {
01486       coarsenElements[i].len=coarsenElements[i/2].len;
01487       coarsenElements[i].elID=coarsenElements[i/2].elID;  // swap bubble with number above  
01488       i/=2;                     
01489     }
01490     coarsenElements[i].elID=eIdx;
01491     coarsenElements[i].len=len; 
01492   }
01493   else {
01494     i = ++refineHeapSize; 
01495     while ((refineElements[i/2].len>=len) && (i != 1))
01496     {
01497       refineElements[i].len=refineElements[i/2].len;
01498       refineElements[i].elID=refineElements[i/2].elID;  // swap bubble with number above  
01499       i/=2;                     
01500     }
01501     refineElements[i].elID=eIdx;
01502     refineElements[i].len=len; 
01503   }
01504 }
01505 
01506 // removes and returns the minimum element of the heap 
01507 int chunk::Delete_Min(int cflag)
01508 {
01509   int Child, i, Min_ID; 
01510 
01511   if (cflag) {
01512     Min_ID=coarsenElements[1].elID;
01513     for (i=1; i*2 <= coarsenHeapSize-1; i=Child)
01514     {
01515       // Find smaller child
01516       Child = i*2;       // child is left child  
01517       if (Child !=coarsenHeapSize)  // If this is not == Size, then right child exists
01518     if (coarsenElements[Child+1].len < coarsenElements[Child].len)
01519       Child++; 
01520 
01521       // Percolate one level
01522       if (coarsenElements[coarsenHeapSize].len >= coarsenElements[Child].len){  
01523     coarsenElements[i].elID = coarsenElements[Child].elID;   
01524     coarsenElements[i].len = coarsenElements[Child].len;
01525       }
01526       else
01527          break; 
01528     }
01529     coarsenElements[i].elID = coarsenElements[coarsenHeapSize].elID;   // Place holding value into the bubble
01530     coarsenElements[i].len = coarsenElements[coarsenHeapSize].len; 
01531     coarsenHeapSize--;
01532     return Min_ID; 
01533   }
01534   else {
01535     Min_ID=refineElements[1].elID;
01536     for (i=1; i*2 <= refineHeapSize-1; i=Child)
01537     {
01538       // Find smaller child
01539       Child = i*2;       // child is left child  
01540       if (Child !=refineHeapSize)  // If this is not == Size, then right child exists
01541     if (refineElements[Child+1].len < refineElements[Child].len)
01542       Child++; 
01543 
01544       // Percolate one level
01545       if (refineElements[refineHeapSize].len >= refineElements[Child].len){  
01546     refineElements[i].elID = refineElements[Child].elID;   
01547     refineElements[i].len = refineElements[Child].len;
01548       }
01549       else
01550          break; 
01551     }
01552     refineElements[i].elID = refineElements[refineHeapSize].elID;   // Place holding value into the bubble
01553     refineElements[i].len = refineElements[refineHeapSize].len; 
01554     refineHeapSize--;
01555     return Min_ID; 
01556   }
01557 }
01558 
01559 
01560 
01561 intMsg *chunk::getBoundary(int edgeIdx)
01562 {
01563   intMsg *im = new intMsg;
01564   im->anInt = theEdges[edgeIdx].getBoundary();
01565   return im;
01566 }
01567 
01568 
01569 void chunk::validCheck()
01570 {
01571   int i;
01572   FEM_Entity *fem_node = meshPtr->lookup(FEM_NODE,"validCheck");
01573   FEM_DataAttribute *validNode = (FEM_DataAttribute *)fem_node->lookup(FEM_ATTRIB_TAG_MAX-1,"validCheck");
01574   AllocTable2d<int> &validNodeTable = validNode->getInt();
01575   /*
01576   for (i=0; i<elementSlots; i++) {
01577     //    CkAssert(theElements[i].isPresent()==);
01578   }
01579   */
01580   for (i=0; i<nodeSlots; i++) {
01581     CkAssert(theNodes[i].isPresent()==validNodeTable[i][0]);
01582   }
01583 }
01584 
01585 intMsg *chunk::neighboring(int idx, elemRef e)
01586 {
01587   intMsg *m= new intMsg;
01588   m->anInt = theElements[idx].neighboring(e);
01589   return m;
01590 }
01591 
01592 intMsg *chunk::safeToCoarsen(int idx, edgeRef e)
01593 {
01594   intMsg *m= new intMsg;
01595   m->anInt = theElements[idx].safeToCoarsen(e);
01596   return m;
01597 }
01598 
01599 #include "refine.def.h"
01600
libs/ck-libs/TMRC2D/tri.C
