Charm++: Charm Source Code Documentation

00001 /* File: fem_adapt.C
00002  * Authors: Nilesh Choudhury, Terry Wilmarth
00003  *
00004  */
00005 
00006 #include "fem_adapt.h"  
00007 
00008 // ======================  BEGIN edge_flip  =================================
00009 /* Perform a Delaunay flip of the edge (n1, n2) returning 1 if
00010    successful, 0 if not (likely due to the edge being on a boundary).
00011    The convexity of the quadrilateral formed by two faces incident to
00012    edge (n1, n2) is assumed. n1 and n2 are assumed to be local to this
00013    chunk.  An adjacency test is performed on n1 and n2 by searching
00014    for an element with edge [n1,n2]. */
00015 int FEM_Adapt::edge_flip(int n1, int n2) 
00016 {
00017   int e1, e2, e1_n1, e1_n2, e1_n3, n3, edge1, e1nbr;
00018   findAdjData(n1, n2, &e1, &e2, &e1_n1, &e1_n2, &e1_n3, &n3, &edge1, &e1nbr);
00019   if (e2 == -1) return 0; // there is no neighbor; can't flip; 
00020   return edge_flip_help(e1, e2, n1, n2, e1_n1, e1_n2, e1_n3, n3, edge1, e1nbr);
00021 }
00022 int FEM_Adapt::edge_flip_help(int e1, int e2, int n1, int n2, int e1_n1, 
00023                               int e1_n2, int e1_n3, int n3, int edge1, 
00024                               int e1nbr) 
00025 {
00026   int e2_n1 = find_local_node_index(e2, n1);
00027   int e2_n2 = find_local_node_index(e2, n2);
00028   int e2_n3 = 3 - e2_n1 - e2_n2;
00029   int mod_edge2 = get_edge_index(e2_n1, e2_n3);
00030   int e2nbr = theMesh->e2e_getNbr(e2, mod_edge2);
00031   int n4 = theMesh->e2n_getNode(e2, e2_n3);
00032   
00033   // Element-to-node updates
00034   theMesh->e2n_setIndex(e1, e1_n2, n4);
00035   theMesh->e2n_setIndex(e2, e2_n1, n3);
00036   // Element-to-element updates
00037   theMesh->e2e_replace(e1nbr, e1, e2);
00038   theMesh->e2e_replace(e2nbr, e2, e1);
00039   theMesh->e2e_replace(e1, e2, e2nbr);
00040   theMesh->e2e_replace(e2, e1, e1nbr);
00041   theMesh->e2e_setIndex(e1, edge1, e2);
00042   theMesh->e2e_setIndex(e2, mod_edge2, e1);
00043   // Node-to-node updates
00044   theMesh->n2n_remove(n1, n2);
00045   theMesh->n2n_remove(n2, n1);
00046   theMesh->n2n_add(n3, n4);
00047   theMesh->n2n_add(n4, n3);
00048   // Node-to-element updates
00049   theMesh->n2e_remove(n1, e2);
00050   theMesh->n2e_remove(n2, e1);
00051   theMesh->n2e_add(n3, e2);
00052   theMesh->n2e_add(n4, e1);
00053   return 1;
00054 }
00055 // ======================  END edge_flip  ===================================
00056 
00057 
00058 // ======================  BEGIN edge_bisect  ===============================
00059 /* Given edge e:(n1, n2), determing the two elements e1:(n1,n2,n3) and 
00060    e2:(n2,n1,n4) adjacent to e, and bisect e, e1 and e2 by adding a node 
00061    n5 to edge e. e1 becomes (n1,n5,n3) and we add e3:(n5,n2,n3), e2 becomes 
00062    (n5,n1,n4) and we add e4:(n2,n5,n4); returns new node n5.
00063 
00064        n3                 n3
00065         o                  o
00066        / \                /|\
00067       /   \              / | \ e1nbr
00068      / e1  \            /e1|e3\
00069     /       \          /   |n5 \
00070 n1 o---------o n2  n1 o----o----o n2
00071     \       /          \   |   / 
00072      \ e2  /            \e2|e4/
00073       \   /              \ | / e2nbr
00074        \ /                \|/
00075         o                  o
00076        n4                 n4
00077 */
00078 int FEM_Adapt::edge_bisect(int n1, int n2) 
00079 {
00080   int e1, e2, e1_n1, e1_n2, e1_n3, n3, edge1, e1nbr;
00081   findAdjData(n1, n2, &e1, &e2, &e1_n1, &e1_n2, &e1_n3, &n3, &edge1, &e1nbr);
00082   return edge_bisect_help(e1, e2, n1, n2, e1_n1, e1_n2, e1_n3, n3,edge1,e1nbr);
00083 }
00084 int FEM_Adapt::edge_bisect_help(int e1, int e2, int n1, int n2, int e1_n1, 
00085                                 int e1_n2, int e1_n3, int n3, int edge1, 
00086                                 int e1nbr)
00087 {
00088   int n5 = newNode();
00089   int e3 = newElement();
00090   theMesh->setMeshSizing(e3, getMeshSizing(e1));
00091   localEdgeBisect(n1, n2, e1, e2, e3, e1_n1, e1_n2, e1_n3, e1nbr, n3, n5);
00092 
00093   if (e2 != -1) {
00094     int e2_n1, e2_n2, e2_n3, n4, edge2, e2nbr;
00095     findAdjData(n1, n2, e2, &e2_n1, &e2_n2, &e2_n3, &n4, &edge2, &e2nbr);
00096     int e4 = newElement();
00097     theMesh->setMeshSizing(e4, getMeshSizing(e2));
00098     localEdgeBisect(n1, n2, e2, e1, e4, e2_n1, e2_n2, e2_n3, e2nbr, n4, n5);
00099     theMesh->n2n_remove(n5, n1);
00100     theMesh->n2n_remove(n5, n2);
00101     theMesh->e2e_replace(e3, -1, e4);
00102     theMesh->e2e_replace(e4, -1, e3);
00103   }
00104   return n5;
00105 }
00106 void FEM_Adapt::localEdgeBisect(int n1, int n2, int e1, int e2, int e3, 
00107                                 int e1n1, int e1n2, int e1n3, int e1nbr, 
00108                                 int n3, int n5)
00109 {
00110   theMesh->e2n_setIndex(e1, e1n2, n5);
00111   theMesh->e2n_setIndex(e3, e1n1, n5);
00112   theMesh->e2n_setIndex(e3, e1n2, n2);
00113   theMesh->e2n_setIndex(e3, e1n3, n3);
00114   theMesh->e2e_replace(e1, e1n3, e3);
00115   theMesh->e2e_replace(e1nbr, e1, e3); // check for -1 and ghost
00116   int nl[3];
00117   nl[get_edge_index(e1n1, e1n2)] = -1;
00118   nl[get_edge_index(e1n2, e1n3)] = e1nbr;
00119   nl[get_edge_index(e1n3, e1n1)] = e1;
00120   theMesh->e2e_setAll(e3, nl);
00121   // update e3 with e4 later
00122   theMesh->n2n_replace(n1, n2, n5);
00123   theMesh->n2n_replace(n2, n1, n5);
00124   theMesh->n2n_add(n3, n5);
00125   theMesh->n2n_add(n5, n1);
00126   theMesh->n2n_add(n5, n2);
00127   theMesh->n2n_add(n5, n3);
00128   // update n5 with n4 later
00129   theMesh->n2e_replace(n2, e1, e3);
00130   theMesh->n2e_add(n3, e3);
00131   theMesh->n2e_add(n5, e1);
00132   theMesh->n2e_add(n5, e3);
00133   // update n5 with e2 & e4 later
00134 }
00135 // ======================  END edge_bisect  ================================
00136 
00137 
00138 // ======================  BEGIN vertex_remove  ============================
00139 /* Inverse of edge bisect, this removes a degree 4 vertex n1 and 2 of its
00140    adjacent elements.  n2 indicates that the two elements removed are
00141    adjacent to edge [n1,n2]. This could be performed with edge_contraction,
00142    but this is a simpler operation. */
00143 int FEM_Adapt::vertex_remove(int n1, int n2)
00144 {
00145   int e1, e2, e1_n1, e1_n2, e1_n3, n3, edge1, e1nbr;
00146   findAdjData(n1, n2, &e1, &e2, &e1_n1, &e1_n2, &e1_n3, &n3, &edge1, &e1nbr);
00147   return vertex_remove_help(e1, e2, n1, n2, e1_n1, e1_n2,e1_n3,n3,edge1,e1nbr);
00148 }
00149 int FEM_Adapt::vertex_remove_help(int e1, int e2, int n1, int n2, int e1_n1,
00150                                   int e1_n2, int e1_n3, int n3, int edge1, 
00151                                   int e1nbr)
00152 {
00153   int e3 = theMesh->e2e_getNbr(e1, get_edge_index(e1_n1, e1_n3));
00154   int e4 = e1nbr;
00155   int e5=-1, e6=-1, n4=-1, n5=-1;
00156   if (e2 > -1) {
00157     int e2_n1, e2_n2, e2_n3;
00158     e2_n1 = find_local_node_index(e2, n1);
00159     e2_n2 = find_local_node_index(e2, n2);
00160     e2_n3 = 3 - e2_n1 - e2_n2;
00161     e5 = theMesh->e2e_getNbr(e2, get_edge_index(e2_n1, e2_n3));
00162     e6 = theMesh->e2e_getNbr(e2, get_edge_index(e2_n2, e2_n3));
00163     n4 = theMesh->e2n_getNode(e2, e2_n3);
00164   }
00165   // find n5
00166   int *nbrNodes, nnsize;
00167   theMesh->n2n_getAll(n1, &nbrNodes, &nnsize);
00168   for (int i=0; i<nnsize; i++) {
00169     if ((nbrNodes[i] != n2) && (nbrNodes[i] != n3) && (nbrNodes[i] != n4)) {
00170       n5 = nbrNodes[i];
00171       break;
00172     }
00173   }
00174   
00175   // Element-to-node updates
00176   theMesh->e2n_replace(e3, n1, n2);
00177   theMesh->e2n_replace(e5, n1, n2);
00178   // Element-to-element updates
00179   theMesh->e2e_replace(e3, e1, e4);
00180   theMesh->e2e_replace(e5, e2, e6);
00181   theMesh->e2e_replace(e4, e1, e3);
00182   theMesh->e2e_replace(e6, e2, e5);
00183   // Node-to-node updates
00184   theMesh->n2n_remove(n3, n1);
00185   theMesh->n2n_remove(n4, n1);
00186   theMesh->n2n_replace(n2, n1, n5);
00187   theMesh->n2n_replace(n5, n1, n2);
00188   // Node-to-element updates
00189   theMesh->n2e_replace(n2, e1, e3);
00190   theMesh->n2e_replace(n2, e2, e5);
00191   theMesh->n2e_remove(n3, e1);
00192   theMesh->n2e_remove(n4, e2);
00193   
00194   deleteNode(n1);
00195   deleteElement(e1);
00196   deleteElement(e2);
00197   return 1;
00198 }
00199 // ======================  END vertex_remove  ==============================
00200   
00201 // edge_contraction and helpers
00202 int FEM_Adapt::edge_contraction(int n1, int n2) 
00203 {
00204   int e1 = theMesh->getElementOnEdge(n1, n2);
00205   if (e1 < 0) {
00206     CkPrintf("ERROR: edge_contraction: no element with edge [%d,%d]\n", n1,n2);
00207     return 0;        
00208   }
00209   return edge_contraction(e1, n1, n2);
00210 }
00211 int FEM_Adapt::edge_contraction(int e1, int n1, int n2) 
00212 {
00213   int e1_n1 = find_local_node_index(e1, n1);
00214   int e1_n2 = find_local_node_index(e1, n2);
00215   int shared_edge = get_edge_index(e1_n1, e1_n2);
00216   int e2 = theMesh->e2e_getNbr(e1, shared_edge); 
00217   return edge_contraction_help(e1, e2, n1, n2, e1_n1, e1_n2);
00218 }
00219 int FEM_Adapt::edge_contraction_help(int e1, int e2, int n1, int n2, int e1_n1,
00220                                      int e1_n2)
00221 {
00222   int e1_n3 = 3 - e1_n1 - e1_n2;
00223   int mod_edge1 = get_edge_index(e1_n1, e1_n3);
00224   int e1nbr1 = theMesh->e2e_getNbr(e1, mod_edge1);
00225   int mod_edge2 = get_edge_index(e1_n2, e1_n3);
00226   int e1nbr2 = theMesh->e2e_getNbr(e1, mod_edge2);
00227   int n3 = theMesh->e2n_getNode(e1, e1_n3);
00228   int e2_n1;
00229   int e2_n2;
00230   int e2_n3;
00231   int mod_edge3;
00232   int e2nbr1;
00233   int mod_edge4;
00234   int e2nbr2;
00235   int n4;
00236   if (e2 > -1) {
00237     e2_n1 = find_local_node_index(e2, n1);
00238     e2_n2 = find_local_node_index(e2, n2);
00239     e2_n3 = 3 - e2_n1 - e2_n2;
00240     mod_edge3 = get_edge_index(e2_n1, e2_n3);
00241     e2nbr1 = theMesh->e2e_getNbr(e2, mod_edge3);
00242     mod_edge4 = get_edge_index(e2_n2, e2_n3);
00243     e2nbr2 = theMesh->e2e_getNbr(e2, mod_edge4);
00244     n4 = theMesh->e2n_getNode(e2, e2_n3);
00245   }
00246 
00247   int *n2_nbrNodes, *n2_nbrElems;
00248   int nnsize, nesize;
00249   theMesh->n2n_getAll(n2, &n2_nbrNodes, &nnsize);
00250   theMesh->n2e_getAll(n2, &n2_nbrElems, &nesize);
00251   
00252   // Element-to-node updates
00253   for (int i=0; i<nesize; i++) {
00254     theMesh->e2n_replace(n2_nbrElems[i], n2, n1);
00255     theMesh->n2e_add(n1, n2_nbrElems[i]);
00256   }
00257   theMesh->n2e_remove(n1, e1); // assume duplicated by loop add above
00258   theMesh->n2e_remove(n1, e2); // assume duplicated by loop add above
00259   // Element-to-element updates
00260   theMesh->e2e_replace(e1nbr1, e1, e1nbr2);
00261   theMesh->e2e_replace(e1nbr2, e1, e1nbr1);
00262   theMesh->e2e_replace(e2nbr1, e2, e2nbr2);
00263   theMesh->e2e_replace(e2nbr2, e2, e2nbr1);
00264   // Node-to-node updates
00265   for (int i=0; i<nnsize; i++) {
00266     if (n2_nbrNodes[i] != n1) {
00267       theMesh->n2n_remove(n2_nbrNodes[i], n1);
00268       theMesh->n2n_replace(n2_nbrNodes[i], n2, n1);
00269       theMesh->n2n_remove(n1, n2_nbrNodes[i]);
00270       theMesh->n2n_add(n1, n2_nbrNodes[i]);
00271     }
00272   }
00273   theMesh->n2n_remove(n1, n2);
00274   theMesh->n2n_remove(n3, n2);
00275   theMesh->n2n_remove(n4, n2);
00276   // Node-to-element updates
00277   theMesh->n2e_remove(n1, e1);
00278   theMesh->n2e_remove(n1, e2);
00279   theMesh->n2e_remove(n3, e1);
00280   theMesh->n2e_remove(n4, e2);
00281   
00282   deleteNode(n2);
00283   deleteElement(e1);
00284   deleteElement(e2);
00285   return 1;
00286 }
00287 
00288 // vertex_split and helpers
00289 int FEM_Adapt::vertex_split(int n, int n1, int n2) 
00290 {
00291   int e1 = theMesh->getElementOnEdge(n, n1);
00292   if (e1 < 0) {
00293     CkPrintf("ERROR: vertex_split: no element with edge [%d,%d]\n", n, n1);
00294     return -1;       
00295   }
00296   int e3 = theMesh->getElementOnEdge(n, n2);
00297   if (e3 < 0) {
00298     CkPrintf("ERROR: vertex_split: no element with edge [%d,%d]\n", n, n2);
00299     return -1;       
00300   }
00301   return vertex_split(n, n1, n2, e1, e3);
00302 }
00303 int FEM_Adapt::vertex_split(int n, int n1, int n2, int e1, int e3)
00304 {
00305   int n_e1 = find_local_node_index(e1, n);
00306   int n1_e1 = find_local_node_index(e1, n1);
00307   int e2 = theMesh->e2e_getNbr(e1, get_edge_index(n_e1, n1_e1));
00308   int n_e3 = find_local_node_index(e3, n);
00309   int n2_e3 = find_local_node_index(e3, n2);
00310   int e4 = theMesh->e2e_getNbr(e3, get_edge_index(n_e3, n2_e3));
00311   if (!check_orientation(e1, e3, n, n1, n2)) {
00312     int tmp = e3;
00313     e3 = e4;
00314     e4 = tmp;
00315     n_e3 = find_local_node_index(e3, n);
00316     n2_e3 = find_local_node_index(e3, n2);
00317   }
00318   int np = newNode();
00319   int e5 = newElement();
00320   int e6 = newElement();
00321   int nnCount=0, neCount=0;
00322   int np_nodes[50], np_elems[50]; // I certainly hope the mesh is not this bad
00323   adj_traverse(n, n1, n2, e2, e4, &nnCount, &neCount, np_nodes, np_elems);
00324 
00325   // Element-to-node updates
00326   int nl[3];
00327   if ((n_e1 < n1_e1) || ((n_e1 == 2) && (n1_e1 == 0))) {
00328     nl[0] = n1; nl[1] = n; nl[2] = np;
00329     theMesh->e2n_setAll(e5, nl);
00330     nl[0] = e1; nl[1] = e6; nl[2] = e2;
00331     theMesh->e2e_setAll(e5, nl);
00332   }
00333   else {
00334     nl[0] = n; nl[1] = n1; nl[2] = np;
00335     theMesh->e2n_setAll(e5, nl);
00336     nl[0] = e1; nl[1] = e2; nl[2] = e6;
00337     theMesh->e2e_setAll(e5, nl);
00338   }
00339   if ((n_e3 < n2_e3) || ((n_e3 == 2) && (n2_e3 == 0))) {
00340     nl[0] = n2; nl[1] = n; nl[2] = np;
00341     theMesh->e2n_setAll(e6, nl);
00342     nl[0] = e3; nl[1] = e5; nl[2] = e4;
00343     theMesh->e2e_setAll(e6, nl);
00344   }
00345   else {
00346     nl[0] = n; nl[1] = n2; nl[2] = np;
00347     theMesh->e2n_setAll(e6, nl);
00348     nl[0] = e3; nl[1] = e4; nl[2] = e5;
00349     theMesh->e2e_setAll(e6, nl);
00350   }
00351   theMesh->e2n_replace(e2, n, np);
00352   theMesh->e2n_replace(e4, n, np);
00353   // Element-to-element updates
00354   theMesh->e2e_replace(e1, e2, e5);
00355   theMesh->e2e_replace(e2, e1, e5);
00356   theMesh->e2e_replace(e3, e4, e6);
00357   theMesh->e2e_replace(e4, e3, e6);
00358   // Node-to-node updates
00359   int i;
00360   for (i=0; i<nnCount; i++) {
00361     printf("np_nodes[%d] = %d\n", i, np_nodes[i]);
00362     theMesh->n2n_remove(n, np_nodes[i]);
00363     theMesh->n2n_remove(np_nodes[i], n);
00364     theMesh->n2n_add(np, np_nodes[i]);
00365     theMesh->n2n_add(np_nodes[i], np);
00366   }
00367   theMesh->n2n_add(n, np);
00368   theMesh->n2n_add(np, n);
00369   theMesh->n2n_add(n, n1);
00370   theMesh->n2n_add(n1, n);
00371   theMesh->n2n_add(n, n2);
00372   theMesh->n2n_add(n2, n);
00373   // Node-to-element updates
00374   for (i=0; i<neCount; i++) { 
00375     theMesh->n2e_remove(n, np_elems[i]);
00376     theMesh->e2n_replace(np_elems[i], n, np);
00377     theMesh->n2e_add(np, np_elems[i]);
00378   }
00379   theMesh->n2e_add(n, e5);
00380   theMesh->n2e_add(n, e6);
00381   theMesh->n2e_add(n1, e5);
00382   theMesh->n2e_add(n2, e6);
00383   theMesh->n2e_add(np, e5);
00384   theMesh->n2e_add(np, e6);
00385   return np;
00386 }
00387 
00388 // element_bisect and helpers
00389 /*
00390 void FEM_Adapt::element_bisect(int e1) 
00391 {
00392   int mn = theMesh->getMarkedNode(e1);
00393   int mnIdx = theMesh->getMarkedNodeIdx(e1);
00394   int n1Idx = (mnIdx+1)%3;
00395   int n1 = theMesh->e2n_getNode(e1, n1Idx);
00396   int n2Idx = (n1Idx+1)%3;
00397   int n2 = theMesh->e2n_getNode(e1, n2Idx);
00398   int e2 = theMesh->e2e_getNbr(e1, get_edge_index(n1Idx, n2Idx));
00399   int mn2 = theMesh->getMarkedNode(e2);
00400   while (!((mn2 != n1) && (mn2 != n2))) {
00401     element_bisect(e2);
00402     mn = theMesh->getMarkedNode(e1);
00403     mnIdx = theMesh->getMarkedNodeIdx(e1);
00404     n1Idx = (mnIdx+1)%3;
00405     n1 = theMesh->e2n_getNode(e1, n1Idx);
00406     n2Idx = (n1Idx+1)%3;
00407     n2 = theMesh->e2n_getNode(e1, n2Idx);
00408     e2 = theMesh->e2e_getNbr(e1, get_edge_index(n1Idx, n2Idx));
00409     mn2 = theMesh->getMarkedNode(e2);
00410   }
00411   (void) edge_bisect(e1, n1, n2);
00412 }
00413 */
00414 
00415 
00416 int FEM_Adapt::newSlot(FEM_DataAttribute *validAttr){
00417         FEM_Entity *entity = validAttr->getEntity();
00418         AllocTable2d<int> *validData = &validAttr->getInt();
00419         int length = validAttr->getLength();
00420         
00421 /*      printf("valid array before new element length %d\n",length);
00422         printValidArray(validAttr);*/
00423         for(int i=0;i<length;i++){
00424                 if((*validData)[i][0] == 0){
00425                   (*validData)[i][0] = 1;
00426 /*                      printf("valid array after new element at %d \n",i);
00427                         printValidArray(validAttr);*/
00428                   return i;
00429                 }
00430         }
00431         entity->setLength(length+1);
00432         validData = &validAttr->getInt();
00433         (*validData)[length][0] = 1;
00434 /*      printf("valid array after new element current length %d\n",validAttr->getLength());
00435         printValidArray(validAttr);*/
00436         return length;
00437 };
00438 
00439 void FEM_Adapt::invalidateSlot(FEM_DataAttribute *validAttr,int slotNumber){
00440         AllocTable2d<int> *validData = &validAttr->getInt();
00441         (*validData)[slotNumber][0] = 0;
00442 };
00443 
00444 int FEM_Adapt::newNode(){
00445         if(nodeValid){
00446                 return  newSlot(nodeValid);
00447         }else{
00448                 return theMesh->node.size();
00449         }
00450 };
00451 
00452 int FEM_Adapt::newElement(){
00453         if(elemValid){
00454                 return newSlot(elemValid);
00455         }else{
00456                 FEM_Elem *elem = (FEM_Elem *)theMesh->lookup(FEM_ELEM,"newElement");
00457                 return elem->size();
00458         }
00459 };
00460 
00461 void FEM_Adapt::deleteNode(int n){
00462   theMesh->n2e_removeAll(n);
00463   theMesh->n2n_removeAll(n);
00464   invalidateSlot(nodeValid,n);
00465 };
00466 
00467 void FEM_Adapt::deleteElement(int e){
00468   theMesh->e2e_removeAll(e);
00469   invalidateSlot(elemValid,e);
00470 };
00471 
00472 void FEM_Adapt::printValidArray(FEM_DataAttribute *validAttr){
00473   FEM_Entity *entity = validAttr->getEntity();
00474   AllocTable2d<int> *validData = &validAttr->getInt();
00475   int length = entity->getMax();
00476   /*    printf("Valid array---\n");
00477         for(int i=0;i<length;i++){
00478         printf("%d %d\n",i,(*validData)[i][0]);
00479         }*/
00480 };
00481 
00482 
00483 // Helpers
00484 
00485 FEM_DataAttribute *FEM_Adapt::validDataFor(int entityNumber){
00486         FEM_Entity *entity = theMesh->lookup(entityNumber,"validDataFor");
00487         FEM_DataAttribute *validAttribute = (FEM_DataAttribute *)entity->lookup(FEM_IS_VALID,"validDataFor");
00488         return validAttribute;
00489 }
00490 
00491 int FEM_Adapt::get_edge_index(int local_node1, int local_node2) 
00492 {
00493   int sum = local_node1 + local_node2;
00494   CkAssert(local_node1 != local_node2);
00495   if (sum == 1) return 0;
00496   else if (sum == 3) return 1;
00497   else if (sum == 2) return 2;
00498   else {
00499     CkPrintf("ERROR: local node pair is strange: [%d,%d]\n", local_node1,
00500             local_node2);
00501     CkAbort("ERROR: local node pair is strange\n");
00502     return -1;
00503   }
00504 }
00505 
00506 int FEM_Adapt::find_local_node_index(int e, int n) {
00507   int result = theMesh->e2n_getIndex(e, n);
00508   if (result < 0) {
00509     CkPrintf("ERROR: node %d not found on element %d\n", n, e);
00510     CkAbort("ERROR: node not found\n");
00511   }
00512   return result;
00513 }
00514 
00515 int FEM_Adapt::check_orientation(int e1, int e3, int n, int n1, int n2)
00516 {
00517   int e1_n = find_local_node_index(e1, n);
00518   int e1_n1 = find_local_node_index(e1, n1);
00519   int e3_n = find_local_node_index(e3, n);
00520   int e3_n2 = find_local_node_index(e3, n2);
00521   
00522   if (((e1_n1 == (e1_n+1)%3) && (e3_n == (e3_n2+1)%3)) ||
00523       ((e1_n == (e1_n1+1)%3) && (e3_n2 == (e3_n+1)%3)))
00524     return 1;
00525   else return 0;
00526 }
00527 
00528 void FEM_Adapt::adj_traverse(int n, int startNode, int stopNode, int startElem,
00529                              int stopElem, int *nn, int *ne, int *nodeList,
00530                              int *elemList)
00531 {
00532   int elm = startElem, nod = startNode;
00533   int nIdx;
00534   (*nn) = 0;
00535   (*ne) = 0;
00536   if (elm == -1) {
00537     nodeList[*nn] = nod; (*nn)++;
00538   }
00539   while ((elm != stopElem) && (elm > -1)) {
00540     nodeList[*nn] = nod; (*nn)++;
00541     elemList[*ne] = elm; (*ne)++;
00542     nIdx = 3 - find_local_node_index(elm,n) - find_local_node_index(elm,nod);
00543     nod = theMesh->e2n_getNode(elm, nIdx);
00544     elm = theMesh->e2e_getNbr(elm, get_edge_index(find_local_node_index(elm,n),
00545                                                   nIdx));
00546   }
00547   if (elm == stopElem) {
00548     nodeList[*nn] = nod; (*nn)++;
00549     elemList[*ne] = elm; (*ne)++;
00550     nodeList[*nn] = stopNode; (*nn)++;
00551   }
00552   else {
00553     nodeList[*nn] = nod; (*nn)++;
00554     int elm = stopElem;
00555     int nod = stopNode;
00556     while (elm > -1) {
00557       nodeList[*nn] = nod; (*nn)++;
00558       elemList[*ne] = elm; (*ne)++;
00559       nIdx = 3 - find_local_node_index(elm,n) - find_local_node_index(elm,nod);
00560       nod = theMesh->e2n_getNode(elm, nIdx);
00561       elm = theMesh->e2e_getNbr(elm, get_edge_index(find_local_node_index(elm,n), nIdx));
00562     }
00563     nodeList[*nn] = nod; (*nn)++;
00564   }
00565 }
00566 
00567 void FEM_Adapt::findAdjData(int n1, int n2, int *e1, int *e2, int *en1, 
00568                              int *en2, int *en3, int *n3, int *edge, int *nbr)
00569 {
00570   (*e1) = theMesh->getElementOnEdge(n1, n2);
00571   CkAssert((*e1) >= 0); // e1 must exist and be local
00572   (*en1) = find_local_node_index((*e1), n1);
00573   (*en2) = find_local_node_index((*e1), n2);
00574   (*en3) = 3 - (*en1) - (*en2);
00575   (*e2) = theMesh->e2e_getNbr((*e1), get_edge_index((*en1), (*en2)));
00576   (*edge) = get_edge_index((*en2), (*en3));
00577   (*nbr) = theMesh->e2e_getNbr((*e1), (*edge));
00578   (*n3) = theMesh->e2n_getNode((*e1), (*en3));
00579 }
00580 
00581 void FEM_Adapt::findAdjData(int n1, int n2, int e2, int *en1, int *en2, 
00582                             int *en3, int *n4, int *edge, int *nbr)
00583 {
00584   (*en1) = find_local_node_index(e2, n1);
00585   (*en2) = find_local_node_index(e2, n2);
00586   (*en3) = 3 - (*en1) - (*en2);
00587   (*edge) = get_edge_index((*en2), (*en3));
00588   (*nbr) = theMesh->e2e_getNbr(e2, (*edge));
00589   (*n4) = theMesh->e2n_getNode(e2, (*en3));
00590 }
00591 
00592 //init the validData method
00593 
00594
libs/ck-libs/fem/fem_adapt.C
