mesh2dual.c

Go to the documentation of this file.
00001 /*
00002     This file is part of magpar.
00003 
00004     Copyright (C) 2002-2010 Werner Scholz
00005 
00006     www:   http://www.magpar.net/
00007     email: magpar(at)magpar.net
00008 
00009     magpar is free software; you can redistribute it and/or modify
00010     it under the terms of the GNU General Public License as published by
00011     the Free Software Foundation; either version 2 of the License, or
00012     (at your option) any later version.
00013 
00014     magpar is distributed in the hope that it will be useful,
00015     but WITHOUT ANY WARRANTY; without even the implied warranty of
00016     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00017     GNU General Public License for more details.
00018 
00019     You should have received a copy of the GNU General Public License
00020     along with magpar; if not, write to the Free Software
00021     Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
00022 */
00023 
00024 static char const Id[] = "$Id: mesh2dual.c 2962 2010-02-04 19:50:44Z scholz $\n\n";
00025 static char const Td[] = "$Today: " __FILE__ "  " __DATE__ "  "  __TIME__  " $\n\n";
00026 
00027 #include "util.h"
00028 
00029 #ifdef METIS
00030 EXTERN_C_BEGIN
00031 #include "metis.h"
00032 EXTERN_C_END
00033 
00034 #ifdef DEBUG
00035 /* for size of connectivity matrix generated by METIS_MeshToNodal */
00036 #define JASIZE 25
00037 #endif
00038 #endif
00039 
00040 int Mesh2VertEle(int n_ele, int n_vert, int *elevert, int **ia, int **ja)
00041 {
00042   MagparFunctionLogBegin;
00043 
00044   /* This function calculates the vertex-element connectivity.
00045    * For each vertex all elements, which share this vertex are stored.
00046    *
00047    * This is different from METIS's MeshToNodal, which calculates
00048    * the connectivity between nodes (i.e. it connects two vertices [nodes] in
00049    * the graph if they are connected by an edge in the mesh).
00050    */
00051 
00052   int *tca;
00053   ierr = PetscMalloc(n_vert*sizeof(int),&tca);CHKERRQ(ierr);
00054   ierr = PetscMemzero(tca,n_vert*sizeof(int));CHKERRQ(ierr);
00055 
00056   /* count number of elements sharing each node */
00057   for (int i=0;i<n_ele;i++) {
00058     for (int j=0;j<NV;j++) {
00059       tca[elevert[NV*i+j]]++;
00060     }
00061   }
00062 
00063   /* create index array */
00064   int *tia;
00065   ierr = PetscMalloc((n_vert+1)*sizeof(int),&tia);CHKERRQ(ierr);
00066   tia[0]=0;
00067   for (int i=1;i<n_vert+1;i++) {
00068     tia[i]=tia[i-1]+tca[i-1];
00069   }
00070 
00071   PetscPrintf(PETSC_COMM_WORLD,
00072     "Allocating %i MB for %i elements in adjacency matrix\n",
00073     tia[n_vert]*sizeof(int)/1024/1024+1,
00074     tia[n_vert]
00075   );
00076 
00077   int *tja;
00078   ierr = PetscMalloc(tia[n_vert]*sizeof(int),&tja);CHKERRQ(ierr);
00079 
00080   ierr = PetscMemzero(tca,n_vert*sizeof(int));CHKERRQ(ierr);
00081 
00082   /* find all elements, which share a node */
00083   for (int i=0;i<n_ele;i++) {
00084     for (int j=0;j<NV;j++) {
00085       int v;
00086       v=elevert[NV*i+j];
00087       tja[tia[v]+tca[v]]=i;
00088       tca[v]++;
00089     }
00090   }
00091 
00092   /* count total number of connectivities */
00093   int cnt=0;
00094   for (int i=0;i<n_vert;i++) {
00095     cnt+=tca[i];
00096     assert(tca[i]==tia[i+1]-tia[i]);
00097   }
00098   assert(cnt==tia[n_vert]);
00099 
00100 #ifdef DEBUG
00101   PetscPrintf(PETSC_COMM_WORLD,"Dual mesh:\n");
00102   PetscPrintf(PETSC_COMM_WORLD,"   node  |   #  | elements; ia\n");
00103 
00104   int cnt2=0;
00105   for (int i=0;i<n_vert;i++) {
00106     PetscPrintf(PETSC_COMM_WORLD,"%7i %6i ",i,tia[i+1]-tia[i]);
00107     for (int j=tia[i];j<tia[i+1];j++) {
00108       PetscPrintf(PETSC_COMM_WORLD,"%7i ",tja[j]);
00109       cnt2++;
00110     }
00111     PetscPrintf(PETSC_COMM_WORLD,"%6i\n",tia[i]);
00112   }
00113   PetscPrintf(PETSC_COMM_WORLD,"%7i %6i %6i\n",n_vert,0,tia[n_vert]);
00114 
00115   assert(cnt==cnt2);
00116 #endif
00117 
00118   ierr = PetscFree(tca);CHKERRQ(ierr);
00119 
00120   *ia=tia;
00121   *ja=tja;
00122 
00123   MagparFunctionLogReturn(0);
00124 }
00125 
00126 
00127 int Mesh2Nodal(int n_ele, int n_vert, int *elevert, int **ia, int **ja)
00128 {
00129   MagparFunctionLogBegin;
00130 
00131   /* This function calculates the vertex-vertex connectivity.
00132    * For each vertex all vertices, which are connected by an edge, are stored.
00133    *
00134    * This is (for magpar's purposes) a full replacement for METIS's MeshToNodal function.
00135    */
00136 
00137   /* obtain vertex-element connectivity first */
00138   int *veia,*veja;
00139   ierr = Mesh2VertEle(n_ele,n_vert,elevert,&veia,&veja);CHKERRQ(ierr);
00140 
00141   /* allocate memory for counters */
00142   int *tca;
00143   ierr = PetscMalloc(n_vert*sizeof(int),&tca);CHKERRQ(ierr);
00144   ierr = PetscMemzero(tca,n_vert*sizeof(int));CHKERRQ(ierr);
00145 
00146   /* allocate temp. memory for neighbors */
00147 #define MAXVERTVERT 10000
00148   int *tj;
00149   ierr = PetscMalloc(MAXVERTVERT*sizeof(int),&tj);CHKERRQ(ierr);
00150 
00151   /* count number of neighboring vertices */
00152   for (int v=0;v<n_vert;v++) {
00153     for (int i=0;i<MAXVERTVERT;i++) {
00154       tj[i]=-1;
00155     }
00156     for (int i=veia[v];i<veia[v+1];i++) {
00157       int el;
00158       el=veja[i];
00159       for (int j=0;j<NV;j++) {
00160         int vel;
00161         vel=elevert[NV*el+j];
00162 
00163         /* skip self */
00164         if (vel==v) continue;
00165 
00166         int found;
00167         found=0;
00168         for (int k=0;k<tca[v];k++) {
00169           if (vel==tj[k]) {
00170             found=1;
00171             break;
00172           }
00173         }
00174         if (!found) {
00175           /* add new vertex */
00176           tj[tca[v]]=vel;
00177           tca[v]++;
00178           /* check for array overrun */
00179           assert(tca[v]<MAXVERTVERT);
00180         }
00181       }
00182     }
00183   }
00184 
00185   ierr = PetscFree(tj);CHKERRQ(ierr);
00186 
00187   /* create index array */
00188   int *tia;
00189   ierr = PetscMalloc((n_vert+1)*sizeof(int),&tia);CHKERRQ(ierr);
00190   tia[0]=0;
00191   for (int i=1;i<n_vert+1;i++) {
00192     tia[i]=tia[i-1]+tca[i-1];
00193   }
00194 
00195   /* reset counter array */
00196   ierr = PetscMemzero(tca,n_vert*sizeof(int));CHKERRQ(ierr);
00197 
00198   PetscPrintf(PETSC_COMM_WORLD,
00199     "Allocating %i MB for %i elements in adjacency matrix\n",
00200     tia[n_vert]*sizeof(int)/1024/1024+1,
00201     tia[n_vert]
00202   );
00203 
00204   int *tja;
00205   ierr = PetscMalloc(tia[n_vert]*sizeof(int),&tja);CHKERRQ(ierr);
00206   for (int i=0;i<tia[n_vert];i++) {
00207     tja[i]=-1;
00208   }
00209 
00210   /* find all neighboring vertices */
00211   for (int v=0;v<n_vert;v++) {
00212     for (int i=veia[v];i<veia[v+1];i++) {
00213       int el;
00214       el=veja[i];
00215       for (int j=0;j<NV;j++) {
00216         int vel;
00217         vel=elevert[NV*el+j];
00218 
00219         /* skip self */
00220         if (vel==v) continue;
00221 
00222         int found;
00223         found=0;
00224         for (int k=tia[v];k<tia[v+1];k++) {
00225           if (vel==tja[k]) {
00226             found=1;
00227             break;
00228           }
00229         }
00230         if (!found) {
00231           /* add new vertex */
00232           tja[tia[v]+tca[v]]=vel;
00233           tca[v]++;
00234           /* check for array overrun */
00235           assert(tca[v]<=tia[v+1]-tia[v]);
00236         }
00237       }
00238     }
00239     assert(tca[v]==tia[v+1]-tia[v]);
00240   }
00241 
00242   /* count total number of connectivities */
00243   int cnt=0;
00244   for (int i=0;i<n_vert;i++) {
00245     cnt+=tca[i];
00246     assert(tca[i]==tia[i+1]-tia[i]);
00247   }
00248   assert(cnt==tia[n_vert]);
00249 
00250 #ifdef DEBUG
00251   PetscPrintf(PETSC_COMM_WORLD,"Nodal mesh:\n");
00252   PetscPrintf(PETSC_COMM_WORLD,"   node  |   #  | nodes; ia\n");
00253 
00254   int cnt2=0;
00255   for (int i=0;i<n_vert;i++) {
00256     PetscPrintf(PETSC_COMM_WORLD,"%7i %6i ",i,tia[i+1]-tia[i]);
00257     for (int j=tia[i];j<tia[i+1];j++) {
00258       PetscPrintf(PETSC_COMM_WORLD,"%7i ",tja[j]);
00259       cnt2++;
00260     }
00261     PetscPrintf(PETSC_COMM_WORLD,"%6i\n",tia[i]);
00262   }
00263   PetscPrintf(PETSC_COMM_WORLD,"%7i %6i %6i\n",n_vert,0,tia[n_vert]);
00264 
00265   assert(cnt==cnt2);
00266 #endif
00267 
00268   ierr = PetscFree(tca);CHKERRQ(ierr);
00269   ierr = PetscFree(veia);CHKERRQ(ierr);
00270   ierr = PetscFree(veja);CHKERRQ(ierr);
00271 
00272   *ia=tia;
00273   *ja=tja;
00274 
00275 #ifdef METIS
00276 #ifdef DEBUG
00277   /* compare with Metis's METIS_MeshToNodal function */
00278   idxtype *mia,*mja;
00279   ierr = PetscMalloc((n_vert+1)*sizeof(idxtype),&mia);CHKERRQ(ierr);
00280   ierr = PetscMalloc(JASIZE*n_vert*sizeof(idxtype),&mja);CHKERRQ(ierr);
00281   int  zero=0,two=2;
00282   METIS_MeshToNodal(
00283     &n_ele,
00284     &n_vert,
00285     elevert,
00286     &two,
00287     &zero,
00288     mia,
00289     mja
00290   );
00291   /* check for array overrun */
00292   assert(mia[n_vert]<JASIZE*n_vert);
00293 
00294   if (tia[n_vert]!=mia[n_vert]) {
00295     PetscPrintf(PETSC_COMM_WORLD,"ERROR1: %i %i %i %i %i\n",n_vert,tia[n_vert],mia[n_vert]);
00296     SETERRQ(PETSC_ERR_LIB,"Inconsistency with METIS detected! Exiting!\n");
00297   }
00298   for (int v=0;v<n_vert+1;v++) {
00299     if (tia[v]!=mia[v]) {
00300       PetscPrintf(PETSC_COMM_WORLD,"ERROR2: %i %i %i\n",v,tia[v],mia[v]);
00301       SETERRQ(PETSC_ERR_LIB,"Inconsistency with METIS detected! Exiting!\n");
00302     }
00303     for (int i=tia[v];i<tia[n_vert];i++) {
00304       if (tja[i]!=mja[i]) {
00305         PetscPrintf(PETSC_COMM_WORLD,"ERROR3: %i %i %i %i %i\n",v,tia[v],mia[v],tja[i],mja[i]);
00306         SETERRQ(PETSC_ERR_LIB,"Inconsistency with METIS detected! Exiting!\n");
00307       }
00308     }
00309   }
00310 
00311   ierr = PetscFree(mia);CHKERRQ(ierr);
00312   ierr = PetscFree(mja);CHKERRQ(ierr);
00313 #endif
00314 #endif
00315 
00316   MagparFunctionLogReturn(0);
00317 }
00318 
00319 
00320 int Mesh2Dual(int n_ele, int n_vert, int *elevert, int **ia, int **ja)
00321 {
00322   MagparFunctionLogBegin;
00323 
00324   /* This function calculates the element-element connectivity.
00325    * For each element all elements, which share a face, are stored.
00326    *
00327    * This is (for magpar's purposes) a full replacement for METIS's MeshToDual function.
00328    */
00329 
00330   /* obtain vertex-element connectivity first */
00331   int *veia,*veja;
00332   ierr = Mesh2VertEle(n_ele,n_vert,elevert,&veia,&veja);CHKERRQ(ierr);
00333 
00334   /* allocate memory for counters */
00335   int *tca;
00336   ierr = PetscMalloc(n_ele*sizeof(int),&tca);CHKERRQ(ierr);
00337   ierr = PetscMemzero(tca,n_ele*sizeof(int));CHKERRQ(ierr);
00338 
00339   /* allocate temp. memory for neighbors */
00340 #define MAXELEELE NF
00341   int *tj;
00342   ierr = PetscMalloc(MAXELEELE*sizeof(int),&tj);CHKERRQ(ierr);
00343 
00344   /* count number of neighboring elements */
00345   for (int el=0;el<n_ele;el++) {
00346     for (int i=0;i<MAXELEELE;i++) {
00347       tj[i]=-1;
00348     }
00349     for (int j=0;j<NV;j++) {
00350       if (tca[el]==NF) break;
00351       int v;
00352       v=elevert[NV*el+j];
00353       for (int i=veia[v];i<veia[v+1];i++) {
00354         if (tca[el]==NF) break;
00355         if (j==0 && tca[el]==NF-1) break;
00356         int nel;
00357         nel=veja[i];
00358         if (nel==el) continue;
00359         int x;
00360         x=tetnb(elevert+NV*el,elevert+NV*nel,PETSC_NULL);
00361         assert(x<NV);
00362         if (x<3) continue;
00363 
00364         int found;
00365         found=0;
00366 
00367         for (int k=0;k<tca[el];k++) {
00368           if (nel==tj[k]) {
00369             found=1;
00370             break;
00371           }
00372         }
00373         if (!found) {
00374           /* add new neighbor */
00375           tj[tca[el]]=nel;
00376           tca[el]++;
00377           /* check for array overrun */
00378           assert(tca[el]<=MAXELEELE);
00379           if (tca[el]==NF) break;
00380         }
00381       }
00382     }
00383   }
00384 
00385   ierr = PetscFree(tj);CHKERRQ(ierr);
00386 
00387   /* create index array */
00388   int *tia;
00389   ierr = PetscMalloc((n_ele+1)*sizeof(int),&tia);CHKERRQ(ierr);
00390   tia[0]=0;
00391   for (int i=1;i<n_ele+1;i++) {
00392     tia[i]=tia[i-1]+tca[i-1];
00393   }
00394 
00395   /* reset counter array */
00396   ierr = PetscMemzero(tca,n_ele*sizeof(int));CHKERRQ(ierr);
00397 
00398   PetscPrintf(PETSC_COMM_WORLD,
00399     "Allocating %i MB for %i elements in adjacency matrix\n",
00400     tia[n_ele]*sizeof(int)/1024/1024+1,
00401     tia[n_ele]
00402   );
00403 
00404   int *tja;
00405   ierr = PetscMalloc(tia[n_ele]*sizeof(int),&tja);CHKERRQ(ierr);
00406   for (int i=0;i<tia[n_ele];i++) {
00407     tja[i]=-1;
00408   }
00409 
00410   /* store all neighboring elements */
00411   for (int el=0;el<n_ele;el++) {
00412     for (int j=0;j<NV;j++) {
00413       if (tca[el]==NF) break;
00414       int v;
00415       v=elevert[NV*el+j];
00416       for (int i=veia[v];i<veia[v+1];i++) {
00417         if (j==0 && tca[el]==NF-1) break;
00418         if (tca[el]==NF) break;
00419         int nel;
00420         nel=veja[i];
00421         if (nel==el) continue;
00422         int x;
00423         x=tetnb(elevert+NV*el,elevert+NV*nel,PETSC_NULL);
00424         assert(x<NV);
00425         if (x<3) continue;
00426 
00427         int found;
00428         found=0;
00429 
00430         for (int k=0;k<tca[el];k++) {
00431           if (nel==tja[tia[el]+k]) {
00432             found=1;
00433             break;
00434           }
00435         }
00436         if (!found) {
00437           /* add new neighbor */
00438           tja[tia[el]+tca[el]]=nel;
00439           tca[el]++;
00440           /* check for array overrun */
00441           assert(tca[el]<=MAXELEELE);
00442           if (tca[el]==NF) break;
00443         }
00444       }
00445     }
00446     assert(tca[el]==tia[el+1]-tia[el]);
00447   }
00448 
00449   /* count total number of connectivities */
00450   int cnt=0;
00451   for (int i=0;i<n_ele;i++) {
00452     cnt+=tca[i];
00453     assert(tca[i]==tia[i+1]-tia[i]);
00454   }
00455   assert(cnt==tia[n_ele]);
00456 
00457 #ifdef DEBUG
00458   PetscPrintf(PETSC_COMM_WORLD,"Element connectivity:\n");
00459   PetscPrintf(PETSC_COMM_WORLD,"   element  |   #  | neighbors; ia\n");
00460 
00461   int cnt2=0;
00462   for (int i=0;i<n_ele;i++) {
00463     PetscPrintf(PETSC_COMM_WORLD,"%7i %6i ",i,tia[i+1]-tia[i]);
00464     for (int j=tia[i];j<tia[i+1];j++) {
00465       PetscPrintf(PETSC_COMM_WORLD,"%7i ",tja[j]);
00466       cnt2++;
00467     }
00468     PetscPrintf(PETSC_COMM_WORLD,"%6i\n",tia[i]);
00469   }
00470   PetscPrintf(PETSC_COMM_WORLD,"%7i %6i %6i\n",n_ele,0,tia[n_ele]);
00471 
00472   assert(cnt==cnt2);
00473 #endif
00474 
00475   ierr = PetscFree(tca);CHKERRQ(ierr);
00476   ierr = PetscFree(veia);CHKERRQ(ierr);
00477   ierr = PetscFree(veja);CHKERRQ(ierr);
00478 
00479   *ia=tia;
00480   *ja=tja;
00481 
00482 #ifdef METIS
00483 #ifdef DEBUG
00484   /* adjacency matrices */
00485   idxtype *mia,*mja;
00486   ierr = PetscMalloc((n_ele+1)*sizeof(idxtype),&mia);CHKERRQ(ierr);
00487   ierr = PetscMalloc(JASIZE*n_ele*sizeof(idxtype),&mja);CHKERRQ(ierr);
00488   int  zero=0,two=2;
00489   METIS_MeshToDual(
00490     &n_ele,
00491     &n_vert,
00492     elevert,
00493     &two,
00494     &zero,
00495     mia,
00496     mja
00497   );
00498   /* check for array overrun */
00499   assert(mia[n_ele]<NV*n_ele);
00500 
00501   if (tia[n_ele]!=mia[n_ele]) {
00502     PetscPrintf(PETSC_COMM_WORLD,"ERROR1: %i %i %i %i %i\n",n_ele,tia[n_ele],mia[n_ele]);
00503     SETERRQ(PETSC_ERR_LIB,"Inconsitency with METIS detected! Exiting!\n");
00504   }
00505   for (int v=0;v<n_ele+1;v++) {
00506     if (tia[v]!=mia[v]) {
00507       PetscPrintf(PETSC_COMM_WORLD,"ERROR2: %i %i %i\n",v,tia[v],mia[v]);
00508       SETERRQ(PETSC_ERR_LIB,"Inconsitency with METIS detected! Exiting!\n");
00509     }
00510 /* METIS uses different ordering!
00511     for (int i=tia[v];i<tia[n_ele];i++) {
00512       if (tja[i]!=mja[i]) {
00513         PetscPrintf(PETSC_COMM_WORLD,"ERROR3: %i %i %i %i %i\n",v,tia[v],mia[v],tja[i],mja[i]);
00514         SETERRQ(PETSC_ERR_LIB,"Inconsitency with METIS detected! Exiting!\n");
00515       }
00516     }
00517 */
00518   }
00519 
00520   ierr = PetscFree(mia);CHKERRQ(ierr);
00521   ierr = PetscFree(mja);CHKERRQ(ierr);
00522 #endif
00523 #endif
00524 
00525 
00526   MagparFunctionLogReturn(0);
00527 }
00528 
00529 
00530 
00531 int Mesh2Dual2(int n_ele, int n_vert, int *elevert, int **ia, int **ja)
00532 {
00533   MagparFunctionLogBegin;
00534 
00535   /* This function calculates the element-element connectivity.
00536    * For each element all elements, which share a node (!), are stored.
00537    *
00538    * This is different from Mesh2Dual!!!
00539    */
00540 
00541   /* obtain vertex-element connectivity first */
00542   int *veia,*veja;
00543   ierr = Mesh2VertEle(n_ele,n_vert,elevert,&veia,&veja);CHKERRQ(ierr);
00544 
00545   /* allocate memory for counters */
00546   int *tca;
00547   ierr = PetscMalloc(n_ele*sizeof(int),&tca);CHKERRQ(ierr);
00548   ierr = PetscMemzero(tca,n_ele*sizeof(int));CHKERRQ(ierr);
00549 
00550   /* allocate temp. memory for neighbors */
00551 #define MAXELEELE NF
00552   int *tj;
00553   ierr = PetscMalloc(MAXELEELE*sizeof(int),&tj);CHKERRQ(ierr);
00554 
00555   /* count number of neighboring elements */
00556   for (int el=0;el<n_ele;el++) {
00557     for (int i=0;i<MAXELEELE;i++) {
00558       tj[i]=-1;
00559     }
00560     for (int j=0;j<NV;j++) {
00561       int v;
00562       v=elevert[NV*el+j];
00563       for (int i=veia[v];i<veia[v+1];i++) {
00564         int nel;
00565         nel=veja[i];
00566         if (nel==el) continue;
00567 #if 0
00568         /* not necessary since all neighbors will be counted */
00569         int x;
00570         x=tetnb(elevert+NV*el,elevert+NV*nel,PETSC_NULL);
00571         assert(x<NV);
00572         if (x<1) continue;
00573 #endif
00574 
00575         int found;
00576         found=0;
00577 
00578         for (int k=0;k<tca[el];k++) {
00579           if (nel==tj[k]) {
00580             found=1;
00581             break;
00582           }
00583         }
00584         if (!found) {
00585           /* add new neighbor */
00586           tj[tca[el]]=nel;
00587           tca[el]++;
00588           /* check for array overrun */
00589           assert(tca[el]<=MAXELEELE);
00590         }
00591       }
00592     }
00593   }
00594 
00595   ierr = PetscFree(tj);CHKERRQ(ierr);
00596 
00597   /* create index array */
00598   int *tia;
00599   ierr = PetscMalloc((n_ele+1)*sizeof(int),&tia);CHKERRQ(ierr);
00600   tia[0]=0;
00601   for (int i=1;i<n_ele+1;i++) {
00602     tia[i]=tia[i-1]+tca[i-1];
00603   }
00604 
00605   /* reset counter array */
00606   ierr = PetscMemzero(tca,n_ele*sizeof(int));CHKERRQ(ierr);
00607 
00608   PetscPrintf(PETSC_COMM_WORLD,
00609     "Allocating %i MB for %i elements in adjacency matrix\n",
00610     tia[n_ele]*sizeof(int)/1024/1024+1,
00611     tia[n_ele]
00612   );
00613 
00614   int *tja;
00615   ierr = PetscMalloc(tia[n_ele]*sizeof(int),&tja);CHKERRQ(ierr);
00616   for (int i=0;i<tia[n_ele];i++) {
00617     tja[i]=-1;
00618   }
00619 
00620   /* store all neighboring elements */
00621   for (int el=0;el<n_ele;el++) {
00622     for (int j=0;j<NV;j++) {
00623       int v;
00624       v=elevert[NV*el+j];
00625       for (int i=veia[v];i<veia[v+1];i++) {
00626         int nel;
00627         nel=veja[i];
00628         if (nel==el) continue;
00629 #if 0
00630         /* not necessary since all neighbors will be counted */
00631         int x;
00632         x=tetnb(elevert+NV*el,elevert+NV*nel,PETSC_NULL);
00633         assert(x<NV);
00634         if (x<1) continue;
00635 #endif
00636 
00637         int found;
00638         found=0;
00639 
00640         for (int k=0;k<tca[el];k++) {
00641           if (nel==tja[tia[el]+k]) {
00642             found=1;
00643             break;
00644           }
00645         }
00646         if (!found) {
00647           /* add new neighbor */
00648           tja[tia[el]+tca[el]]=nel;
00649           tca[el]++;
00650           /* check for array overrun */
00651           assert(tca[el]<=MAXELEELE);
00652         }
00653       }
00654     }
00655     assert(tca[el]==tia[el+1]-tia[el]);
00656   }
00657 
00658   /* count total number of connectivities */
00659   int cnt=0;
00660   for (int i=0;i<n_ele;i++) {
00661     cnt+=tca[i];
00662     assert(tca[i]==tia[i+1]-tia[i]);
00663   }
00664   assert(cnt==tia[n_ele]);
00665 
00666 #ifdef DEBUG
00667   PetscPrintf(PETSC_COMM_WORLD,"Element connectivity:\n");
00668   PetscPrintf(PETSC_COMM_WORLD,"   element  |   #  | neighbors; ia\n");
00669 
00670   int cnt2=0;
00671   for (int i=0;i<n_ele;i++) {
00672     PetscPrintf(PETSC_COMM_WORLD,"%7i %6i ",i,tia[i+1]-tia[i]);
00673     for (int j=tia[i];j<tia[i+1];j++) {
00674       PetscPrintf(PETSC_COMM_WORLD,"%7i ",tja[j]);
00675       cnt2++;
00676     }
00677     PetscPrintf(PETSC_COMM_WORLD,"%6i\n",tia[i]);
00678   }
00679   PetscPrintf(PETSC_COMM_WORLD,"%7i %6i %6i\n",n_ele,0,tia[n_ele]);
00680 
00681   assert(cnt==cnt2);
00682 #endif
00683 
00684   ierr = PetscFree(tca);CHKERRQ(ierr);
00685   ierr = PetscFree(veia);CHKERRQ(ierr);
00686   ierr = PetscFree(veja);CHKERRQ(ierr);
00687 
00688   *ia=tia;
00689   *ja=tja;
00690 
00691 #ifdef METIS
00692 #ifdef DEBUG
00693   /* adjacency matrices */
00694   idxtype *mia,*mja;
00695   ierr = PetscMalloc((n_ele+1)*sizeof(idxtype),&mia);CHKERRQ(ierr);
00696   ierr = PetscMalloc(JASIZE*n_ele*sizeof(idxtype),&mja);CHKERRQ(ierr);
00697   int  zero=0,two=2;
00698   METIS_MeshToDual(
00699     &n_ele,
00700     &n_vert,
00701     elevert,
00702     &two,
00703     &zero,
00704     mia,
00705     mja
00706   );
00707   /* check for array overrun */
00708   assert(mia[n_ele]<NV*n_ele);
00709 
00710   if (tia[n_ele]!=mia[n_ele]) {
00711     PetscPrintf(PETSC_COMM_WORLD,"ERROR1: %i %i %i %i %i\n",n_ele,tia[n_ele],mia[n_ele]);
00712     SETERRQ(PETSC_ERR_LIB,"Inconsitency with METIS detected! Exiting!\n");
00713   }
00714   for (int v=0;v<n_ele+1;v++) {
00715     if (tia[v]!=mia[v]) {
00716       PetscPrintf(PETSC_COMM_WORLD,"ERROR2: %i %i %i\n",v,tia[v],mia[v]);
00717       SETERRQ(PETSC_ERR_LIB,"Inconsitency with METIS detected! Exiting!\n");
00718     }
00719 /* METIS uses different ordering!
00720     for (int i=tia[v];i<tia[n_ele];i++) {
00721       if (tja[i]!=mja[i]) {
00722         PetscPrintf(PETSC_COMM_WORLD,"ERROR3: %i %i %i %i %i\n",v,tia[v],mia[v],tja[i],mja[i]);
00723         SETERRQ(PETSC_ERR_LIB,"Inconsitency with METIS detected! Exiting!\n");
00724       }
00725     }
00726 */
00727   }
00728 
00729   ierr = PetscFree(mia);CHKERRQ(ierr);
00730   ierr = PetscFree(mja);CHKERRQ(ierr);
00731 #endif
00732 #endif
00733 
00734 
00735   MagparFunctionLogReturn(0);
00736 }
00737 

magpar - Parallel Finite Element Micromagnetics Package
Copyright (C) 2002-2010 Werner Scholz