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   else
00247     {
00248       CkAbort("ERROR: e2 <= -1 case not handled\n");
00249     }
00250   int *n2_nbrNodes, *n2_nbrElems;
00251   int nnsize, nesize;
00252   theMesh->n2n_getAll(n2, &n2_nbrNodes, &nnsize);
00253   theMesh->n2e_getAll(n2, &n2_nbrElems, &nesize);
00254   
00255   // Element-to-node updates
00256   for (int i=0; i<nesize; i++) {
00257     theMesh->e2n_replace(n2_nbrElems[i], n2, n1);
00258     theMesh->n2e_add(n1, n2_nbrElems[i]);
00259   }
00260   theMesh->n2e_remove(n1, e1); // assume duplicated by loop add above
00261   theMesh->n2e_remove(n1, e2); // assume duplicated by loop add above
00262   // Element-to-element updates
00263   theMesh->e2e_replace(e1nbr1, e1, e1nbr2);
00264   theMesh->e2e_replace(e1nbr2, e1, e1nbr1);
00265   theMesh->e2e_replace(e2nbr1, e2, e2nbr2);
00266   theMesh->e2e_replace(e2nbr2, e2, e2nbr1);
00267   // Node-to-node updates
00268   for (int i=0; i<nnsize; i++) {
00269     if (n2_nbrNodes[i] != n1) {
00270       theMesh->n2n_remove(n2_nbrNodes[i], n1);
00271       theMesh->n2n_replace(n2_nbrNodes[i], n2, n1);
00272       theMesh->n2n_remove(n1, n2_nbrNodes[i]);
00273       theMesh->n2n_add(n1, n2_nbrNodes[i]);
00274     }
00275   }
00276   theMesh->n2n_remove(n1, n2);
00277   theMesh->n2n_remove(n3, n2);
00278   theMesh->n2n_remove(n4, n2);
00279   // Node-to-element updates
00280   theMesh->n2e_remove(n1, e1);
00281   theMesh->n2e_remove(n1, e2);
00282   theMesh->n2e_remove(n3, e1);
00283   theMesh->n2e_remove(n4, e2);
00284   
00285   deleteNode(n2);
00286   deleteElement(e1);
00287   deleteElement(e2);
00288   return 1;
00289 }
00290 
00291 // vertex_split and helpers
00292 int FEM_Adapt::vertex_split(int n, int n1, int n2) 
00293 {
00294   int e1 = theMesh->getElementOnEdge(n, n1);
00295   if (e1 < 0) {
00296     CkPrintf("ERROR: vertex_split: no element with edge [%d,%d]\n", n, n1);
00297     return -1;       
00298   }
00299   int e3 = theMesh->getElementOnEdge(n, n2);
00300   if (e3 < 0) {
00301     CkPrintf("ERROR: vertex_split: no element with edge [%d,%d]\n", n, n2);
00302     return -1;       
00303   }
00304   return vertex_split(n, n1, n2, e1, e3);
00305 }
00306 int FEM_Adapt::vertex_split(int n, int n1, int n2, int e1, int e3)
00307 {
00308   int n_e1 = find_local_node_index(e1, n);
00309   int n1_e1 = find_local_node_index(e1, n1);
00310   int e2 = theMesh->e2e_getNbr(e1, get_edge_index(n_e1, n1_e1));
00311   int n_e3 = find_local_node_index(e3, n);
00312   int n2_e3 = find_local_node_index(e3, n2);
00313   int e4 = theMesh->e2e_getNbr(e3, get_edge_index(n_e3, n2_e3));
00314   if (!check_orientation(e1, e3, n, n1, n2)) {
00315     int tmp = e3;
00316     e3 = e4;
00317     e4 = tmp;
00318     n_e3 = find_local_node_index(e3, n);
00319     n2_e3 = find_local_node_index(e3, n2);
00320   }
00321   int np = newNode();
00322   int e5 = newElement();
00323   int e6 = newElement();
00324   int nnCount=0, neCount=0;
00325   int np_nodes[50], np_elems[50]; // I certainly hope the mesh is not this bad
00326   adj_traverse(n, n1, n2, e2, e4, &nnCount, &neCount, np_nodes, np_elems);
00327 
00328   // Element-to-node updates
00329   int nl[3];
00330   if ((n_e1 < n1_e1) || ((n_e1 == 2) && (n1_e1 == 0))) {
00331     nl[0] = n1; nl[1] = n; nl[2] = np;
00332     theMesh->e2n_setAll(e5, nl);
00333     nl[0] = e1; nl[1] = e6; nl[2] = e2;
00334     theMesh->e2e_setAll(e5, nl);
00335   }
00336   else {
00337     nl[0] = n; nl[1] = n1; nl[2] = np;
00338     theMesh->e2n_setAll(e5, nl);
00339     nl[0] = e1; nl[1] = e2; nl[2] = e6;
00340     theMesh->e2e_setAll(e5, nl);
00341   }
00342   if ((n_e3 < n2_e3) || ((n_e3 == 2) && (n2_e3 == 0))) {
00343     nl[0] = n2; nl[1] = n; nl[2] = np;
00344     theMesh->e2n_setAll(e6, nl);
00345     nl[0] = e3; nl[1] = e5; nl[2] = e4;
00346     theMesh->e2e_setAll(e6, nl);
00347   }
00348   else {
00349     nl[0] = n; nl[1] = n2; nl[2] = np;
00350     theMesh->e2n_setAll(e6, nl);
00351     nl[0] = e3; nl[1] = e4; nl[2] = e5;
00352     theMesh->e2e_setAll(e6, nl);
00353   }
00354   theMesh->e2n_replace(e2, n, np);
00355   theMesh->e2n_replace(e4, n, np);
00356   // Element-to-element updates
00357   theMesh->e2e_replace(e1, e2, e5);
00358   theMesh->e2e_replace(e2, e1, e5);
00359   theMesh->e2e_replace(e3, e4, e6);
00360   theMesh->e2e_replace(e4, e3, e6);
00361   // Node-to-node updates
00362   int i;
00363   for (i=0; i<nnCount; i++) {
00364     printf("np_nodes[%d] = %d\n", i, np_nodes[i]);
00365     theMesh->n2n_remove(n, np_nodes[i]);
00366     theMesh->n2n_remove(np_nodes[i], n);
00367     theMesh->n2n_add(np, np_nodes[i]);
00368     theMesh->n2n_add(np_nodes[i], np);
00369   }
00370   theMesh->n2n_add(n, np);
00371   theMesh->n2n_add(np, n);
00372   theMesh->n2n_add(n, n1);
00373   theMesh->n2n_add(n1, n);
00374   theMesh->n2n_add(n, n2);
00375   theMesh->n2n_add(n2, n);
00376   // Node-to-element updates
00377   for (i=0; i<neCount; i++) { 
00378     theMesh->n2e_remove(n, np_elems[i]);
00379     theMesh->e2n_replace(np_elems[i], n, np);
00380     theMesh->n2e_add(np, np_elems[i]);
00381   }
00382   theMesh->n2e_add(n, e5);
00383   theMesh->n2e_add(n, e6);
00384   theMesh->n2e_add(n1, e5);
00385   theMesh->n2e_add(n2, e6);
00386   theMesh->n2e_add(np, e5);
00387   theMesh->n2e_add(np, e6);
00388   return np;
00389 }
00390 
00391 // element_bisect and helpers
00392 /*
00393 void FEM_Adapt::element_bisect(int e1) 
00394 {
00395   int mn = theMesh->getMarkedNode(e1);
00396   int mnIdx = theMesh->getMarkedNodeIdx(e1);
00397   int n1Idx = (mnIdx+1)%3;
00398   int n1 = theMesh->e2n_getNode(e1, n1Idx);
00399   int n2Idx = (n1Idx+1)%3;
00400   int n2 = theMesh->e2n_getNode(e1, n2Idx);
00401   int e2 = theMesh->e2e_getNbr(e1, get_edge_index(n1Idx, n2Idx));
00402   int mn2 = theMesh->getMarkedNode(e2);
00403   while (!((mn2 != n1) && (mn2 != n2))) {
00404     element_bisect(e2);
00405     mn = theMesh->getMarkedNode(e1);
00406     mnIdx = theMesh->getMarkedNodeIdx(e1);
00407     n1Idx = (mnIdx+1)%3;
00408     n1 = theMesh->e2n_getNode(e1, n1Idx);
00409     n2Idx = (n1Idx+1)%3;
00410     n2 = theMesh->e2n_getNode(e1, n2Idx);
00411     e2 = theMesh->e2e_getNbr(e1, get_edge_index(n1Idx, n2Idx));
00412     mn2 = theMesh->getMarkedNode(e2);
00413   }
00414   (void) edge_bisect(e1, n1, n2);
00415 }
00416 */
00417 
00418 
00419 int FEM_Adapt::newSlot(FEM_DataAttribute *validAttr){
00420     FEM_Entity *entity = validAttr->getEntity();
00421     AllocTable2d<int> *validData = &validAttr->getInt();
00422     int length = validAttr->getLength();
00423     
00424 /*  printf("valid array before new element length %d\n",length);
00425     printValidArray(validAttr);*/
00426     for(int i=0;i<length;i++){
00427         if((*validData)[i][0] == 0){
00428           (*validData)[i][0] = 1;
00429 /*          printf("valid array after new element at %d \n",i);
00430             printValidArray(validAttr);*/
00431           return i;
00432         }
00433     }
00434     entity->setLength(length+1);
00435     validData = &validAttr->getInt();
00436     (*validData)[length][0] = 1;
00437 /*  printf("valid array after new element current length %d\n",validAttr->getLength());
00438     printValidArray(validAttr);*/
00439     return length;
00440 };
00441 
00442 void FEM_Adapt::invalidateSlot(FEM_DataAttribute *validAttr,int slotNumber){
00443     AllocTable2d<int> *validData = &validAttr->getInt();
00444     (*validData)[slotNumber][0] = 0;
00445 };
00446 
00447 int FEM_Adapt::newNode(){
00448     if(nodeValid){
00449         return  newSlot(nodeValid);
00450     }else{
00451         return theMesh->node.size();
00452     }
00453 };
00454 
00455 int FEM_Adapt::newElement(){
00456     if(elemValid){
00457         return newSlot(elemValid);
00458     }else{
00459         FEM_Elem *elem = (FEM_Elem *)theMesh->lookup(FEM_ELEM,"newElement");
00460         return elem->size();
00461     }
00462 };
00463 
00464 void FEM_Adapt::deleteNode(int n){
00465   theMesh->n2e_removeAll(n);
00466   theMesh->n2n_removeAll(n);
00467   invalidateSlot(nodeValid,n);
00468 };
00469 
00470 void FEM_Adapt::deleteElement(int e){
00471   theMesh->e2e_removeAll(e);
00472   invalidateSlot(elemValid,e);
00473 };
00474 
00475 void FEM_Adapt::printValidArray(FEM_DataAttribute *validAttr){
00476   FEM_Entity *entity = validAttr->getEntity();
00477   AllocTable2d<int> *validData = &validAttr->getInt();
00478   int length = entity->getMax();
00479   /*    printf("Valid array---\n");
00480     for(int i=0;i<length;i++){
00481     printf("%d %d\n",i,(*validData)[i][0]);
00482     }*/
00483 };
00484 
00485 
00486 // Helpers
00487 
00488 FEM_DataAttribute *FEM_Adapt::validDataFor(int entityNumber){
00489     FEM_Entity *entity = theMesh->lookup(entityNumber,"validDataFor");
00490     FEM_DataAttribute *validAttribute = (FEM_DataAttribute *)entity->lookup(FEM_IS_VALID,"validDataFor");
00491     return validAttribute;
00492 }
00493 
00494 int FEM_Adapt::get_edge_index(int local_node1, int local_node2) 
00495 {
00496   int sum = local_node1 + local_node2;
00497   CkAssert(local_node1 != local_node2);
00498   if (sum == 1) return 0;
00499   else if (sum == 3) return 1;
00500   else if (sum == 2) return 2;
00501   else {
00502     CkPrintf("ERROR: local node pair is strange: [%d,%d]\n", local_node1,
00503         local_node2);
00504     CkAbort("ERROR: local node pair is strange\n");
00505     return -1;
00506   }
00507 }
00508 
00509 int FEM_Adapt::find_local_node_index(int e, int n) {
00510   int result = theMesh->e2n_getIndex(e, n);
00511   if (result < 0) {
00512     CkPrintf("ERROR: node %d not found on element %d\n", n, e);
00513     CkAbort("ERROR: node not found\n");
00514   }
00515   return result;
00516 }
00517 
00518 int FEM_Adapt::check_orientation(int e1, int e3, int n, int n1, int n2)
00519 {
00520   int e1_n = find_local_node_index(e1, n);
00521   int e1_n1 = find_local_node_index(e1, n1);
00522   int e3_n = find_local_node_index(e3, n);
00523   int e3_n2 = find_local_node_index(e3, n2);
00524   
00525   if (((e1_n1 == (e1_n+1)%3) && (e3_n == (e3_n2+1)%3)) ||
00526       ((e1_n == (e1_n1+1)%3) && (e3_n2 == (e3_n+1)%3)))
00527     return 1;
00528   else return 0;
00529 }
00530 
00531 void FEM_Adapt::adj_traverse(int n, int startNode, int stopNode, int startElem,
00532                  int stopElem, int *nn, int *ne, int *nodeList,
00533                  int *elemList)
00534 {
00535   int elm = startElem, nod = startNode;
00536   int nIdx;
00537   (*nn) = 0;
00538   (*ne) = 0;
00539   if (elm == -1) {
00540     nodeList[*nn] = nod; (*nn)++;
00541   }
00542   while ((elm != stopElem) && (elm > -1)) {
00543     nodeList[*nn] = nod; (*nn)++;
00544     elemList[*ne] = elm; (*ne)++;
00545     nIdx = 3 - find_local_node_index(elm,n) - find_local_node_index(elm,nod);
00546     nod = theMesh->e2n_getNode(elm, nIdx);
00547     elm = theMesh->e2e_getNbr(elm, get_edge_index(find_local_node_index(elm,n),
00548                           nIdx));
00549   }
00550   if (elm == stopElem) {
00551     nodeList[*nn] = nod; (*nn)++;
00552     elemList[*ne] = elm; (*ne)++;
00553     nodeList[*nn] = stopNode; (*nn)++;
00554   }
00555   else {
00556     nodeList[*nn] = nod; (*nn)++;
00557     int elm = stopElem;
00558     int nod = stopNode;
00559     while (elm > -1) {
00560       nodeList[*nn] = nod; (*nn)++;
00561       elemList[*ne] = elm; (*ne)++;
00562       nIdx = 3 - find_local_node_index(elm,n) - find_local_node_index(elm,nod);
00563       nod = theMesh->e2n_getNode(elm, nIdx);
00564       elm = theMesh->e2e_getNbr(elm, get_edge_index(find_local_node_index(elm,n), nIdx));
00565     }
00566     nodeList[*nn] = nod; (*nn)++;
00567   }
00568 }
00569 
00570 void FEM_Adapt::findAdjData(int n1, int n2, int *e1, int *e2, int *en1, 
00571                  int *en2, int *en3, int *n3, int *edge, int *nbr)
00572 {
00573   (*e1) = theMesh->getElementOnEdge(n1, n2);
00574   CkAssert((*e1) >= 0); // e1 must exist and be local
00575   (*en1) = find_local_node_index((*e1), n1);
00576   (*en2) = find_local_node_index((*e1), n2);
00577   (*en3) = 3 - (*en1) - (*en2);
00578   (*e2) = theMesh->e2e_getNbr((*e1), get_edge_index((*en1), (*en2)));
00579   (*edge) = get_edge_index((*en2), (*en3));
00580   (*nbr) = theMesh->e2e_getNbr((*e1), (*edge));
00581   (*n3) = theMesh->e2n_getNode((*e1), (*en3));
00582 }
00583 
00584 void FEM_Adapt::findAdjData(int n1, int n2, int e2, int *en1, int *en2, 
00585                 int *en3, int *n4, int *edge, int *nbr)
00586 {
00587   (*en1) = find_local_node_index(e2, n1);
00588   (*en2) = find_local_node_index(e2, n2);
00589   (*en3) = 3 - (*en1) - (*en2);
00590   (*edge) = get_edge_index((*en2), (*en3));
00591   (*nbr) = theMesh->e2e_getNbr(e2, (*edge));
00592   (*n4) = theMesh->e2n_getNode(e2, (*en3));
00593 }
00594 
00595 //init the validData method
00596 
00597
libs/ck-libs/fem/fem_adapt.C
