43557f60bbb829075492cc960fff7044b8a1c308
[charm.git] / src / libs / ck-libs / ParFUM / adapt.C
1 /* File: fem_adapt_new.C
2  * Authors: Nilesh Choudhury, Terry Wilmarth
3  *
4  */
5
6 #include "ParFUM.h"
7 #include "ParFUM_internals.h"
8
9 //#define DEBUG_1
10
11 /** Makes use of ordering of nodes in e1 to check is e3 is on the same side
12     of the path of edges (n1, n) and (n, n2) */
13 int FEM_Adapt::check_orientation(int e1, int e3, int n, int n1, int n2)
14 {
15   int e1_n = find_local_node_index(e1, n);
16   int e1_n1 = find_local_node_index(e1, n1);
17   int e3_n = find_local_node_index(e3, n);
18   int e3_n2 = find_local_node_index(e3, n2);
19   
20   if (((e1_n1 == (e1_n+1)%3) && (e3_n == (e3_n2+1)%3)) ||
21       ((e1_n == (e1_n1+1)%3) && (e3_n2 == (e3_n+1)%3)))
22     return 1;
23   else return 0;
24 }
25
26 /** Given two element-local node numberings (i.e. 0, 1, 2 for triangular 
27     elements), calculate an element-local edge numbering (also 0, 1, or 2
28     for triangular elements)
29 */
30 int FEM_Adapt::get_edge_index(int local_node1, int local_node2) 
31 {
32   int sum = local_node1 + local_node2;
33   CkAssert(local_node1 != local_node2);
34   if (sum == 1) return 0;
35   else if (sum == 3) return 1;
36   else if (sum == 2) return 2;
37   else {
38     CkPrintf("ERROR: local node pair is strange: [%d,%d]\n", local_node1,
39             local_node2);
40     CkAbort("ERROR: local node pair is strange\n");
41     return -1;
42   }
43 }
44
45 /** Given a chunk-local element number e and a chunk-local node number n,
46     determine the element-local node numbering for node n on element e **/
47 int FEM_Adapt::find_local_node_index(int e, int n) {
48   int result = theMesh->e2n_getIndex(e, n);
49   if (result < 0) {
50     CkPrintf("ERROR: node %d not found on element %d\n", n, e);
51     CkAbort("ERROR: node not found\n");
52   }
53   return result;
54 }
55
56
57
58 /** Extract elements adjacent to edge [n1,n2] along with element-local node
59     numberings and nodes opposite input edge **/
60 int FEM_Adapt::findAdjData(int n1, int n2, int *e1, int *e2, int *e1n1, 
61                             int *e1n2, int *e1n3, int *e2n1, int *e2n2, 
62                             int *e2n3, int *n3, int *n4)
63 {
64   // Set some default values in case e1 is not there
65   (*e1n1) = (*e1n2) = (*e1n3) = (*n3) = -1;
66
67   //if n1,n2 is not an edge return 
68   if(n1<0 || n2<0) return -1;
69   if(theMesh->node.is_valid(n1)==0 || theMesh->node.is_valid(n2)==0) return -1;
70   if(theMesh->n2n_exists(n1,n2)!=1 || theMesh->n2n_exists(n2,n1)!=1) return -1; 
71
72   (*e1) = theMesh->getElementOnEdge(n1, n2); // assumed to return local element
73   if ((*e1) == -1) {
74     CkPrintf("[%d]Warning: No Element on edge %d->%d\n",theMod->idx,n1,n2);
75     return -1;
76   }
77   (*e1n1) = find_local_node_index((*e1), n1);
78   (*e1n2) = find_local_node_index((*e1), n2);
79   (*e1n3) = 3 - (*e1n1) - (*e1n2);
80   (*n3) = theMesh->e2n_getNode((*e1), (*e1n3));
81   (*e2) = theMesh->e2e_getNbr((*e1), get_edge_index((*e1n1), (*e1n2)));
82   // Set some default values in case e2 is not there
83   (*e2n1) = (*e2n2) = (*e2n3) = (*n4) = -1;
84   if ((*e2) != -1) { // e2 exists
85     (*e2n1) = find_local_node_index((*e2), n1);
86     (*e2n2) = find_local_node_index((*e2), n2);
87     (*e2n3) = 3 - (*e2n1) - (*e2n2);
88     //if ((*e2) > -1) { // if e2 is a ghost, there is no e2n data
89     (*n4) = theMesh->e2n_getNode((*e2), (*e2n3));
90     //}
91   }
92   if(*n3 == *n4) {
93     CkPrintf("[%d]Warning: Identical elements %d:(%d,%d,%d) & %d:(%d,%d,%d)\n",theMod->idx,*e1,n1,n2,*n3,*e2,n1,n2,*n4);
94     return -1;
95   }
96   return 1;
97 }
98
99 int FEM_Adapt::e2n_getNot(int e, int n1, int n2) {
100   int eConn[3];
101   theMesh->e2n_getAll(e, eConn);
102   for (int i=0; i<3; i++) 
103     if ((eConn[i] != n2) && (eConn[i] != n1)) return eConn[i];
104   return -1; //should never come here
105 }
106
107 int FEM_Adapt::n2e_exists(int n, int e) {
108   int *nConn, nSz;
109   theMesh->n2e_getAll(n, nConn, nSz);
110   for (int i=0; i<nSz; i++) {
111     if (nConn[i] == e) {
112       if(nSz!=0) free(nConn);
113       return 1;
114     }
115   }
116   if(nSz!=0) free(nConn);
117   return 0;
118 }
119
120 int FEM_Adapt::findElementWithNodes(int n1, int n2, int n3) {
121   int *nConn, nSz;
122   int ret = -1;
123   theMesh->n2e_getAll(n1, nConn, nSz);
124   for (int i=0; i<nSz; i++) {
125     if ((n2e_exists(n2, nConn[i])) && (n2e_exists(n3, nConn[i]))) {
126       ret = nConn[i];
127       break;
128     }
129   }
130   if(nSz!=0) free(nConn);
131   return ret; //should never come here
132 }
133
134
135
136 int FEM_Adapt::getSharedNodeIdxl(int n, int chk) {
137   return theMod->getfmUtil()->exists_in_IDXL(theMesh, n, chk, 0, -1);
138 }
139
140 int FEM_Adapt::getGhostNodeIdxl(int n, int chk) { 
141   return theMod->getfmUtil()->exists_in_IDXL(theMesh, n, chk, 2, -1);
142 }
143
144 int FEM_Adapt::getGhostElementIdxl(int e, int chk) { 
145   return theMod->getfmUtil()->exists_in_IDXL(theMesh, e, chk, 4, 0);
146 }
147
148
149
150 void FEM_Adapt::printAdjacencies(int *nodes, int numNodes, int *elems, int numElems) {
151
152   for(int i=0; i<numNodes; i++) {
153     if(nodes[i] == -1) continue;
154     theMod->getfmUtil()->FEM_Print_n2e(theMesh, nodes[i]);
155     theMod->getfmUtil()->FEM_Print_n2n(theMesh, nodes[i]);
156   }
157   for(int i=0; i<numElems; i++) {
158     if(elems[i] == -1) continue;
159     theMod->getfmUtil()->FEM_Print_e2n(theMesh, elems[i]);
160     theMod->getfmUtil()->FEM_Print_e2e(theMesh, elems[i]);
161   }
162
163   return;
164 }
165
166
167
168 bool FEM_Adapt::isFixedNode(int n1) {
169   for(int i=0; i<theMod->fmfixedNodes.size(); i++) {
170     if(theMod->fmfixedNodes[i]==n1) return true;
171   }
172   return false;
173 }
174
175 /** A node is a corner if it is connected to two different boundaries and it is 
176     on the boundary
177 */
178 bool FEM_Adapt::isCorner(int n1) {
179   //if it has at least two adjacent nodes on different boundaries and the edges are boundaries
180   int *n1AdjNodes;
181   int n1NumNodes=0;
182   int n1_bound, n2_bound;
183   FEM_Mesh_dataP(theMesh, FEM_NODE, FEM_BOUNDARY, &n1_bound, n1, 1 , FEM_INT, 1);
184   if(n1_bound==0) return false; //it is internal
185   theMesh->n2n_getAll(n1, n1AdjNodes, n1NumNodes);
186   for (int i=0; i<n1NumNodes; i++) {
187     int n2 = n1AdjNodes[i];
188     if(FEM_Is_ghost_index(n2)) {
189       int numchunks;
190       IDXL_Share **chunks1;
191       theMod->fmUtil->getChunkNos(0,n2,&numchunks,&chunks1);
192       int index = theMod->idx;
193       CkAssert(numchunks>0);
194       int chk = chunks1[0]->chk;
195       int ghostidx = theMod->fmUtil->exists_in_IDXL(theMesh,n2,chk,2);
196       intMsg *im = meshMod[chk].getRemoteBound(index,ghostidx);
197       n2_bound = im->i;
198       for(int j=0; j<numchunks; j++) {
199         delete chunks1[j];
200       }
201       delete im;
202       if(numchunks>0) free(chunks1);
203     }
204     else {
205       FEM_Mesh_dataP(theMesh, FEM_NODE, FEM_BOUNDARY, &n2_bound, n2, 1 , FEM_INT, 1);
206     }
207     if(n2_bound == 0) continue;
208     if(n1_bound != n2_bound) {
209       if(isEdgeBoundary(n1,n2) && abs(n1_bound)>abs(n2_bound)) {
210         if(n1NumNodes!=0) delete[] n1AdjNodes;
211         return true;
212       }
213     }
214   }
215   if(n1NumNodes!=0) delete[] n1AdjNodes;
216   return false;
217 }
218
219 bool FEM_Adapt::isEdgeBoundary(int n1, int n2) {
220   int *n1AdjElems, *n2AdjElems;
221   int n1NumElems=0, n2NumElems=0;
222   int ret = 0;
223   //find the number of elements this edge belongs to
224   theMesh->n2e_getAll(n1, n1AdjElems, n1NumElems);
225   theMesh->n2e_getAll(n2, n2AdjElems, n2NumElems);
226   for(int k=0; k<n1NumElems; k++) {
227     for (int j=0; j<n2NumElems; j++) {
228       if (n1AdjElems[k] == n2AdjElems[j]) {
229         if(n1AdjElems[k] != -1) {
230           ret++;
231         }
232       }
233     }
234   }
235   if(n1NumElems!=0) delete[] n1AdjElems;
236   if(n2NumElems!=0) delete[] n2AdjElems;
237   if(ret==1) return true;
238   return false;
239 }
240
241
242
243 // ======================  BEGIN edge_flip  =================================
244 /** Perform a Delaunay flip of the edge (n1, n2) returning 1 if successful, 0 if
245    not (likely due to the edge being on a boundary). The convexity of the 
246    quadrilateral formed by two faces incident to edge (n1, n2) is assumed. n1 
247    and n2 are assumed to be local to this chunk.  An adjacency test is 
248    performed on n1 and n2 by searching for an element with edge [n1,n2].
249
250        n3                 n3
251         o                  o
252        / \                /|\
253       /   \              / | \
254      /     \            /  |  \
255     /       \          /   |   \
256 n1 o---------o n2  n1 o    |    o n2
257     \       /          \   |   / 
258      \     /            \  |  /
259       \   /              \ | /
260        \ /                \|/
261         o                  o
262        n4                 n4
263 */
264 int FEM_Adapt::edge_flip_help(int e1, int e2, int n1, int n2, int e1_n1, 
265                               int e1_n2, int e1_n3, int n3, int n4, int *locknodes) 
266 {
267   int numNodes = 4;
268   int numElems = 2;
269   int lockelems[2];
270   int elemConn[3];
271   locknodes[0] = n1;
272   locknodes[1] = n2;
273   locknodes[2] = n3;
274   locknodes[3] = n4;
275   lockelems[0] = e1;
276   lockelems[1] = e2;
277   if(n1 < 0 || n2 < 0) {
278     return -1;
279   }
280   int index = theMod->getIdx();
281   bool flag = theMod->fmAdaptAlgs->controlQualityF(n1,n2,n3,n4);
282   if(flag) return -1;
283   int e1Topurge = e1;
284   int e2Topurge = e2;
285
286 #ifdef DEBUG_1
287   CkPrintf("Flipping edge %d->%d on chunk %d\n", n1, n2, theMod->getfmUtil()->getIdx());
288 #endif
289   //FEM_Modify_Lock(theMesh, locknodes, numNodes, lockelems, numElems);
290   //if any of the two elements is remote, eat those
291   if(n3 < 0) {
292     e1Topurge = theMod->fmUtil->eatIntoElement(e1);
293     theMesh->e2n_getAll(e1Topurge,elemConn);
294     for(int i=0; i<3; i++) {
295       if(elemConn[i]!=n1 && elemConn[i]!=n2) {
296         n3 = elemConn[i];
297       }
298     }
299     locknodes[2] = n3;
300   }
301   if(n4 < 0) {
302     e2Topurge = theMod->fmUtil->eatIntoElement(e2);
303     theMesh->e2n_getAll(e2Topurge,elemConn);
304     for(int i=0; i<3; i++) {
305       if(elemConn[i]!=n1 && elemConn[i]!=n2) {
306         n4 = elemConn[i];
307       }
308     }
309     locknodes[3] = n4;
310   }
311
312   FEM_remove_element(theMesh,e1Topurge,0,0);
313   FEM_remove_element(theMesh,e2Topurge,0,0);
314   // add n1, n3, n4
315   elemConn[e1_n1] = n1;  elemConn[e1_n2] = n4;  elemConn[e1_n3] = n3;
316   lockelems[0] = FEM_add_element(theMesh, elemConn, 3, 0, index);
317   //the attributes should really be interpolated, i.e. on both new elems,
318   //the values should be an average of the previous two elements
319   theMod->fmUtil->copyElemData(0,e1Topurge,lockelems[0]);
320   // add n2, n3, n4
321   elemConn[e1_n1] = n4;  elemConn[e1_n2] = n2;  elemConn[e1_n3] = n3;
322   lockelems[1] = FEM_add_element(theMesh, elemConn, 3, 0, index);
323   theMod->fmUtil->copyElemData(0,e1Topurge,lockelems[1]); //both of the new elements copy from one element
324   //purge the two elements
325   FEM_purge_element(theMesh,e1Topurge,0);
326   FEM_purge_element(theMesh,e2Topurge,0);
327
328   //get rid of some unnecessary ghost node sends
329   for(int i=0; i<4;i++) {
330     int nodeToUpdate = -1;
331     if(i==0) nodeToUpdate = n1;
332     else if(i==1) nodeToUpdate = n2;
333     else if(i==2) nodeToUpdate = n3;
334     else if(i==3) nodeToUpdate = n4;
335     //if any of the chunks sharing this node sends this as a ghost, then all of them have to
336     //so find out the set of chunks I need to send this as a ghost to
337     //collect info from each of the shared chunks, do a union of all these chunks
338     //send this updated list to everyone.
339     //if anyone needs to add or delete some ghosts, they will
340     int *chkl, numchkl=0;
341     CkVec<int> finalchkl;
342     theMod->fmUtil->findGhostSend(nodeToUpdate, chkl, numchkl);
343     for(int j=0; j<numchkl; j++) {
344       finalchkl.push_back(chkl[j]);
345     }
346     if(numchkl>0) delete[] chkl;
347
348     const IDXL_Rec *irec=theMesh->node.shared.getRec(nodeToUpdate);
349     int numchunks=0;
350     int *chunks1, *inds1;
351     if(irec) {
352       numchunks = irec->getShared();
353       chunks1 = new int[numchunks];
354       inds1 = new int[numchunks];
355       for(int j=0; j<numchunks; j++) {
356         chunks1[j] = irec->getChk(j);
357         inds1[j] = irec->getIdx(j);
358       }
359     }
360     for(int j=0; j<numchunks; j++) {
361       findgsMsg *fmsg = meshMod[chunks1[j]].findghostsend(index,inds1[j]);
362       if(fmsg->numchks>0) {
363         for(int k=0; k<fmsg->numchks; k++) {
364           bool shouldbeadded = true;
365           for(int l=0; l<finalchkl.size(); l++) {
366             if(fmsg->chunks[k]==finalchkl[l]) {
367               shouldbeadded = false;
368               break;
369             }
370           }
371           if(shouldbeadded) finalchkl.push_back(fmsg->chunks[k]);
372         }
373       }
374       delete fmsg;
375     }
376
377     int *finall, numfinall=finalchkl.size();
378     if(numfinall>0) finall = new int[numfinall];
379     for(int j=0; j<numfinall; j++) finall[j] = finalchkl[j];
380     finalchkl.free();
381
382     theMod->fmUtil->UpdateGhostSend(nodeToUpdate, finall, numfinall);
383     for(int j=0; j<numchunks; j++) {
384       verifyghostsendMsg *vmsg = new(numfinall)verifyghostsendMsg();
385       vmsg->fromChk = index;
386       vmsg->sharedIdx = inds1[j];
387       vmsg->numchks = numfinall;
388       for(int k=0; k<numfinall; k++) vmsg->chunks[k] = finall[k];
389       meshMod[chunks1[j]].updateghostsend(vmsg);
390     }
391     if(numfinall>0) delete[] finall;
392     if(numchunks>0) {
393       delete[] chunks1;
394       delete[] inds1;
395     }
396   }
397   //make sure that it always comes here, don't return with unlocking
398   return 1; //return newNode;
399 }
400 // ======================  END edge_flip  ===================================
401
402
403 // ======================  BEGIN edge_bisect  ===============================
404 /** Given edge e:(n1, n2), remove the two elements (n1,n2,n3) and 
405    (n2,n1,n4) adjacent to e, and bisect e by adding node 
406    n5. Add elements (n1,n5,n3), (n5,n2,n3), (n5,n1,n4) and (n2,n5,n4); 
407    returns new node n5.
408
409        n3                 n3
410         o                  o
411        / \                /|\
412       /   \              / | \
413      /     \            /  |  \
414     /       \          /   |n5 \
415 n1 o---------o n2  n1 o----o----o n2
416     \       /          \   |   / 
417      \     /            \  |  /
418       \   /              \ | /
419        \ /                \|/
420         o                  o
421        n4                 n4
422 */
423 int FEM_Adapt::edge_bisect_help(int e1, int e2, int n1, int n2, int e1_n1, 
424                                 int e1_n2, int e1_n3, int e2_n1, int e2_n2, 
425                                 int e2_n3, int n3, int n4)
426 {
427   int n5;
428   int numNodes = 4;
429   int numElems = 2;
430   int numNodesNew = 5;
431   int numElemsNew = 4;
432   int locknodes[5];
433   int lockelems[4];
434   int elemConn[3];
435
436   locknodes[0] = n1;
437   locknodes[1] = n2;
438   locknodes[2] = n3;
439   locknodes[3] = n4;
440   locknodes[4] = -1;
441   lockelems[0] = e1;
442   lockelems[1] = e2;
443   lockelems[2] = -1;
444   lockelems[3] = -1;
445
446   //FEM_Modify_Lock(theMesh, locknodes, numNodes, lockelems, numElems);
447
448   int e1chunk=-1, e2chunk=-1, e3chunk=-1, e4chunk=-1, n5chunk=-1;
449   int index = theMod->getIdx();
450
451 #ifdef DEBUG_1
452   CkPrintf("Bisect edge %d->%d on chunk %d\n", n1, n2, theMod->getfmUtil()->getIdx());
453 #endif
454   //verify quality
455   bool flag = theMod->fmAdaptAlgs->controlQualityR(n1,n2,n3,n4);
456   if(flag) return -1;
457
458   //add node
459   if(e1==-1) e1chunk=-1;
460   else if(e1>=0) e1chunk=index;
461   else {
462     int ghostid = FEM_To_ghost_index(e1);
463     const IDXL_Rec *irec = theMesh->elem[0].ghost->ghostRecv.getRec(ghostid);
464     CkAssert(irec->getShared()==1);
465     e1chunk = irec->getChk(0);
466   }
467   if(e2==-1) e2chunk=-1;
468   else if(e2>=0) e2chunk=index;
469   else {
470     int ghostid = FEM_To_ghost_index(e2);
471     const IDXL_Rec *irec = theMesh->elem[0].ghost->ghostRecv.getRec(ghostid);
472     CkAssert(irec->getShared()==1);
473     e2chunk = irec->getChk(0);
474   }
475   int adjnodes[2];
476   adjnodes[0] = n1;
477   adjnodes[1] = n2;
478   int *chunks;
479   int numChunks=0;
480   int forceshared = 0;
481   if(e1chunk==e2chunk || (e1chunk==-1 || e2chunk==-1)) {
482     forceshared = -1;
483     numChunks = 1;
484     chunks = new int[1];
485     if(e1chunk!=-1) chunks[0] = e1chunk;
486     else chunks[0] = e2chunk;
487   }
488   else {
489     numChunks = 2;
490     chunks = new int[2];
491     chunks[0] = e1chunk;
492     chunks[1] = e2chunk;
493   }
494   n5 = FEM_add_node(theMesh,adjnodes,2,chunks,numChunks,forceshared);
495   delete[] chunks;
496   //lock this node immediately
497   FEM_Modify_LockN(theMesh, n5, 0);
498
499   //remove elements
500   e1chunk = FEM_remove_element(theMesh, e1, 0); 
501   e2chunk = FEM_remove_element(theMesh, e2, 0);  // assumes intelligent behavior when no e2 exists
502
503   // hmm... if e2 is a ghost and we remove it and create all the new elements
504   // locally, then we don't really need to add a *shared* node
505   //but we are not moving chunk boundaries for bisect
506   if(e1chunk==-1 || e2chunk==-1) {
507     //it is fine, let it continue
508     e4chunk = e2chunk;
509     e3chunk = e1chunk;
510   }
511   else if(e1chunk==e2chunk && e1chunk!=index) {
512     n5chunk = e1chunk;
513     e4chunk = e2chunk;
514     e3chunk = e1chunk;
515   }
516   else {
517     //there can be a lot of conditions, but we do not have to do aything special now
518     n5chunk = -1;
519     e4chunk = e2chunk;
520     e3chunk = e1chunk;
521   }
522
523   // add n1, n5, n3
524   elemConn[e1_n1] = n1;  elemConn[e1_n2] = n5;  elemConn[e1_n3] = n3;
525   lockelems[0] = FEM_add_element(theMesh, elemConn, 3, 0, e1chunk);
526   theMod->fmUtil->copyElemData(0,e1,lockelems[0]);
527   // add n2, n5, n3
528   elemConn[e1_n1] = n5;  elemConn[e1_n2] = n2;  elemConn[e1_n3] = n3;
529   lockelems[1] = FEM_add_element(theMesh, elemConn, 3, 0, e3chunk);
530   theMod->fmUtil->copyElemData(0,e1,lockelems[1]);
531   if (e2 != -1) { // e2 exists
532     // add n1, n5, n4
533     elemConn[e2_n1] = n1;  elemConn[e2_n2] = n5;  elemConn[e2_n3] = n4;
534     lockelems[2] = FEM_add_element(theMesh, elemConn, 3, 0, e2chunk);
535     theMod->fmUtil->copyElemData(0,e2,lockelems[2]);
536     // add n2, n5, n4
537     elemConn[e2_n1] = n5;  elemConn[e2_n2] = n2;  elemConn[e2_n3] = n4;
538     lockelems[3] = FEM_add_element(theMesh, elemConn, 3, 0, e4chunk);
539     theMod->fmUtil->copyElemData(0,e2,lockelems[3]);
540   }
541
542   FEM_purge_element(theMesh,e1,0);
543   FEM_purge_element(theMesh,e2,0);
544
545   //get rid of some unnecessary ghost node sends
546   for(int i=0; i<4;i++) {
547     int nodeToUpdate = -1;
548     if(i==0) nodeToUpdate = n1;
549     else if(i==1) nodeToUpdate = n2;
550     else if(i==2) nodeToUpdate = n3;
551     else if(i==3) nodeToUpdate = n4;
552     //if any of the chunks sharing this node sends this as a ghost, then all of them have to
553     //so find out the set of chunks I need to send this as a ghost to
554     //collect info from each of the shared chunks, do a union of all these chunks
555     //send this updated list to everyone.
556     //if anyone needs to add or delete some ghosts, they will
557     int *chkl, numchkl=0;
558     CkVec<int> finalchkl;
559     theMod->fmUtil->findGhostSend(nodeToUpdate, chkl, numchkl);
560     for(int j=0; j<numchkl; j++) {
561       finalchkl.push_back(chkl[j]);
562     }
563     if(numchkl>0) delete[] chkl;
564
565     const IDXL_Rec *irec=theMesh->node.shared.getRec(nodeToUpdate);
566     int numchunks=0;
567     int *chunks1, *inds1;
568     if(irec) {
569       numchunks = irec->getShared();
570       chunks1 = new int[numchunks];
571       inds1 = new int[numchunks];
572       for(int j=0; j<numchunks; j++) {
573         chunks1[j] = irec->getChk(j);
574         inds1[j] = irec->getIdx(j);
575       }
576     }
577     for(int j=0; j<numchunks; j++) {
578       findgsMsg *fmsg = meshMod[chunks1[j]].findghostsend(index,inds1[j]);
579       if(fmsg->numchks>0) {
580         for(int k=0; k<fmsg->numchks; k++) {
581           bool shouldbeadded = true;
582           for(int l=0; l<finalchkl.size(); l++) {
583             if(fmsg->chunks[k]==finalchkl[l]) {
584               shouldbeadded = false;
585               break;
586             }
587           }
588           if(shouldbeadded) finalchkl.push_back(fmsg->chunks[k]);
589         }
590       }
591       delete fmsg;
592     }
593
594     int *finall, numfinall=finalchkl.size();
595     if(numfinall>0) finall = new int[numfinall];
596     for(int j=0; j<numfinall; j++) finall[j] = finalchkl[j];
597     finalchkl.free();
598
599     theMod->fmUtil->UpdateGhostSend(nodeToUpdate, finall, numfinall);
600     for(int j=0; j<numchunks; j++) {
601       verifyghostsendMsg *vmsg = new(numfinall)verifyghostsendMsg();
602       vmsg->fromChk = index;
603       vmsg->sharedIdx = inds1[j];
604       vmsg->numchks = numfinall;
605       for(int k=0; k<numfinall; k++) vmsg->chunks[k] = finall[k];
606       meshMod[chunks1[j]].updateghostsend(vmsg);
607     }
608     if(numfinall>0) delete[] finall;
609     if(numchunks>0) {
610       delete[] chunks1;
611       delete[] inds1;
612     }
613   }
614
615   FEM_Modify_UnlockN(theMesh, n5, 0);
616   return n5;
617 }
618 // ======================  END edge_bisect  ================================
619
620
621 // ======================  BEGIN vertex_remove  ============================
622 /** Inverse of edge bisect, this removes a degree 4 vertex n1 and 2 of its
623    adjacent elements.  n2 indicates that the two elements removed are
624    adjacent to edge [n1,n2]. This could be performed with edge_contraction,
625    but this is a simpler operation. 
626
627          n3                  n3        
628           o                   o        
629          /|\                 / \       
630         / | \               /   \      
631        /  |  \             /     \     
632       /   |n1 \           /       \    
633   n5 o----o----o n2   n5 o---------o n2
634       \   |   /           \       /    
635        \  |  /             \     /     
636         \ | /               \   /      
637          \|/                 \ /       
638           o                   o        
639          n4                  n4        
640 */
641 int FEM_Adapt::vertex_remove_help(int e1, int e2, int n1, int n2, int e1_n1, 
642                                   int e1_n2, int e1_n3, int e2_n1, int e2_n2, 
643                                   int e2_n3, int n3, int n4, int n5)
644 {
645   int numNodes = 5;
646   int numElems = 4;
647   int numNodesNew = 4;
648   int numElemsNew = 2;
649   int locknodes[5];
650   int lockelems[4];
651   int elemConn[3];
652
653   locknodes[0] = n2;
654   locknodes[1] = n3;
655   locknodes[2] = n4;
656   locknodes[3] = n5;
657   locknodes[4] = n1;
658   lockelems[0] = e1;
659   lockelems[1] = e2;
660   lockelems[2] = -1;
661   lockelems[3] = -1;
662
663   int e3 = theMesh->e2e_getNbr(e1, get_edge_index(e1_n1, e1_n3));
664   int e4 = -1;
665   lockelems[2] = e3;
666   if (e3 != -1) {
667     if (e2 != -1) {
668       e4 = theMesh->e2e_getNbr(e2, get_edge_index(e2_n1, e2_n3));
669       lockelems[3] = e4;
670       if(e4 == -1 ) {
671         return 0;
672       }
673     }
674     //FEM_Modify_Lock(theMesh, locknodes, numNodes, lockelems, numElems);
675
676     int e1chunk=-1, e2chunk=-1, e3chunk=-1, e4chunk=-1;
677     int index = theMod->getIdx();
678
679 #ifdef DEBUG_1
680     CkPrintf("Vertex Remove edge %d->%d on chunk %d\n", n1, n2, theMod->getfmUtil()->getIdx());
681 #endif
682     
683     e1chunk = FEM_remove_element(theMesh, e1, 0);
684     e3chunk = FEM_remove_element(theMesh, e3, 0);
685     if (e2 != -1) {
686       e2chunk = FEM_remove_element(theMesh, e2, 0);
687       e4chunk = FEM_remove_element(theMesh, e4, 0);
688     }
689     FEM_remove_node(theMesh, n1);
690     
691     // add n2, n5, n3
692     elemConn[e1_n1] = n2;  elemConn[e1_n2] = n3;  elemConn[e1_n3] = n5;
693     lockelems[0] = FEM_add_element(theMesh, elemConn, 3, 0, e1chunk);
694     if (e2 != -1) {
695       // add n2, n5, n4
696       elemConn[e2_n1] = n5;  elemConn[e2_n2] = n4;  elemConn[e2_n3] = n2;
697       lockelems[1] = FEM_add_element(theMesh, elemConn, 3, 0, e2chunk);
698     }
699
700     return 1;
701   }
702
703   return 0;
704 }
705 // ======================  END vertex_remove  ==============================
706   
707 // ======================  BEGIN vertex_split =================================
708 /** Given a node n and two adjacent nodes n1 and n2, split n into two nodes n 
709    and np such that the edges to the neighbors n1 and n2 expand into two new 
710    elements (n, np, n1) and (np, n, n2); return the id of the newly created 
711    node np
712
713     n1              n1             
714      o               o             
715      |              / \            
716      |             /   \           
717    \ | /      \   /     \   /      
718     \|/        \ /       \ /       
719      o n      n o---------o np
720     /|\        / \       / \       
721    / | \      /   \     /   \      
722      |             \   /           
723      |              \ /            
724      o               o             
725     n2              n2             
726 */
727 /* This function has some undefined characteristics.. please do not use it as of now*/
728 int FEM_Adapt::vertex_split(int n, int n1, int n2) 
729 {
730   if ((n < 0) || ((n1 <= -1) && (n2 <= -1)))
731     CkAbort("FEM_Adapt::vertex_split: n and at least one of its neighbor must be local to this chunk; n1 and n2 must both exist\n");
732   int locknodes[2];
733   locknodes[0] = n1; locknodes[1] = n2;
734   FEM_Modify_Lock(theMesh, locknodes, 2, locknodes, 0);
735   int e1 = theMesh->getElementOnEdge(n, n1);
736   if (e1 == -1) {
737     FEM_Modify_Unlock(theMesh);
738     return -1;       
739   }
740   int e3 = theMesh->getElementOnEdge(n, n2);
741   if (e3 == -1) {
742     FEM_Modify_Unlock(theMesh);
743     return -1;       
744   }
745   int ret = vertex_split_help(n, n1, n2, e1, e3);
746   FEM_Modify_Unlock(theMesh);
747   return ret;
748 }
749
750 int FEM_Adapt::vertex_split_help(int n, int n1, int n2, int e1, int e3)
751 {
752   int e1_n = find_local_node_index(e1, n);
753   int e1_n1 = find_local_node_index(e1, n1);
754   int e2 = theMesh->e2e_getNbr(e1, get_edge_index(e1_n, e1_n1));
755   int e3_n = find_local_node_index(e3, n);
756   int e3_n2 = find_local_node_index(e3, n2);
757   int e4 = theMesh->e2e_getNbr(e3, get_edge_index(e3_n, e3_n2));
758   if (!check_orientation(e1, e3, n, n1, n2)) {
759     int tmp = e3;
760     e3 = e4;
761     e4 = tmp;
762     e3_n = find_local_node_index(e3, n);
763     e3_n2 = find_local_node_index(e3, n2);
764   }
765
766   int locknodes[4];
767   locknodes[0] = n1;
768   locknodes[1] = n;
769   locknodes[2] = n2;
770   locknodes[3] = -1;
771   int lockelems[6];
772   lockelems[0] = e1;
773   lockelems[1] = e2;
774   lockelems[2] = e3;
775   lockelems[3] = e4;
776   lockelems[4] = -1;
777   lockelems[5] = -1;
778   int elemConn[3];
779
780   //FEM_Modify_Lock(theMesh, locknodes, 4, lockelems, 6);
781 #ifdef DEBUG_1
782   CkPrintf("Vertex Split, %d-%d-%d on chunk %d\n", n1, n, n2, theMod->getfmUtil()->getIdx());
783 #endif
784   int adjnodes[2];
785   adjnodes[0] = n; //looks like it will never be shared, since according to later code, all n1, n & n2 should be local.. appears to be not correct
786   //the new node will be shared to wahtever the old node was shared to, we'll do this later
787   int *chunks;
788   int numChunks = 0;
789   int np = FEM_add_node(theMesh,adjnodes,1,chunks,numChunks,0);
790   locknodes[3] = np;
791
792   int current, next, nt, nl, eknp, eknt, eknl;
793   // traverse elements on one side of n starting with e2
794   current = e2;
795   nt = n1;
796   while ((current != e3) && (current != -1)) { 
797     eknp = find_local_node_index(current, n);
798     eknt = find_local_node_index(current, nt);
799     eknl = 3 - eknp - eknt;
800     next = theMesh->e2e_getNbr(current, get_edge_index(eknp, eknl));
801     nl = theMesh->e2n_getNode(current, eknl);
802     FEM_remove_element(theMesh, current, 0);
803     // add nl, nt, np
804     elemConn[eknp] = np; elemConn[eknt] = nt; elemConn[eknl] = nl;
805     int newelem = FEM_add_element(theMesh, elemConn, 3, 0);
806     nt = nl;
807     current = next;
808   }
809   if (current == -1) { // didn't make it all the way around
810     // traverse elements on one side of n starting with e4
811     current = e4;
812     nt = n2;
813     while ((current != e1) && (current != -1)) {
814       eknp = find_local_node_index(current, n);
815       eknt = find_local_node_index(current, nt);
816       eknl = 3 - eknp - eknt;
817       next = theMesh->e2e_getNbr(current, get_edge_index(eknp, eknl));
818       nl = theMesh->e2n_getNode(current, eknl);
819       FEM_remove_element(theMesh, current, 0);
820       // add nl, nt, np
821       elemConn[eknp] = np; elemConn[eknt] = nt; elemConn[eknl] = nl;
822       int newelem = FEM_add_element(theMesh, elemConn, 3, 0);
823       nt = nl;
824       current = next;
825     }
826   }
827
828   // add n, n1, np
829   elemConn[e1_n] = n; elemConn[e1_n1] = n1; elemConn[3 - e1_n - e1_n1] = np;
830   lockelems[4] = FEM_add_element(theMesh, elemConn, 3, 0);
831   // add n, n2, np
832   elemConn[e3_n] = n; elemConn[e3_n2] = n2; elemConn[3 - e3_n - e3_n2] = np;
833   lockelems[5] = FEM_add_element(theMesh, elemConn, 3, 0);
834
835   return np;
836 }
837 // ======================  END vertex_split ===================