00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024 static char const Id[] = "$Id: facnb.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 "init.h"
00028 #include "util/util.h"
00029
00030 #define USEQSORT
00031 #ifdef USEQSORT
00032 #include <stdlib.h>
00033
00034 int *facvert;
00035
00036 int CompareFaces(const void* v1,const void* v2) {
00037 const int* t1=(const int*)v1;
00038 const int* t2=(const int*)v2;
00039
00040 if (facvert[NN* *t1+0]<facvert[NN* *t2+0]) return -1;
00041 if (facvert[NN* *t1+0]>facvert[NN* *t2+0]) return 1;
00042
00043 if (facvert[NN* *t1+1]<facvert[NN* *t2+1]) return -1;
00044 if (facvert[NN* *t1+1]>facvert[NN* *t2+1]) return -1;
00045
00046 if (facvert[NN* *t1+2]<facvert[NN* *t2+2]) return -1;
00047 if (facvert[NN* *t1+2]>facvert[NN* *t2+2]) return -1;
00048
00049
00050 assert(0==1);
00051 return 0;
00052 }
00053 #endif
00054
00055 int FacNB(GridData *gdata)
00056 {
00057 int t_facvert[NF*NV];
00058 int t_elevs[NV];
00059
00060 MagparFunctionLogBegin;
00061
00062 int rank,size;
00063 ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRQ(ierr);
00064 ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRQ(ierr);
00065
00066 if (rank) {
00067 gdata->ln_bnd_fac=0;
00068 MagparFunctionLogReturn(0);
00069 }
00070
00071
00072
00073
00074
00075
00076 int *ia,*ja;
00077 ierr = Mesh2Dual(gdata->n_ele,gdata->n_vert,gdata->elevert,&ia,&ja);CHKERRQ(ierr);
00078
00079 for (int i=1;i<gdata->n_ele+1;i++) {
00080 if (ia[i]-ia[i-1]>NF) {
00081 SETERRQ2(PETSC_ERR_MEMC,"Found element with %i neighbors (>NF)!\n",ia[i]-ia[i-1],NF);
00082 }
00083 }
00084 assert(ia[gdata->n_ele]<NF*gdata->n_ele);
00085
00086
00087 assert(ia[gdata->n_ele]<NF*gdata->n_ele);
00088
00089 int *elefac;
00090 ierr = PetscMalloc(gdata->n_ele*NF*sizeof(int),&elefac);CHKERRQ(ierr);
00091
00092 int *facnb;
00093 ierr = PetscMalloc(NF*2*gdata->n_ele*sizeof(int),&facnb);CHKERRQ(ierr);
00094
00095 for (int i=0;i<NF*2*gdata->n_ele;i++) {
00096 facnb[i]=C_UNK;
00097 }
00098
00099
00100
00101
00102
00103
00104
00105
00106
00107 int n_fac;
00108 n_fac=NF*gdata->n_ele-ia[gdata->n_ele]/2;
00109
00110 #ifndef USEQSORT
00111 int *facvert;
00112 #endif
00113 ierr = PetscMalloc(n_fac*NN*sizeof(int),&facvert);CHKERRQ(ierr);
00114
00115 int t_facmax;
00116 t_facmax=0;
00117
00118 ierr = PetscGetTime(&t_t2);CHKERRQ(ierr);
00119
00120 PetscPrintf(PETSC_COMM_WORLD,"Searching for surface faces and nodes...\n");
00121 for (int i=0;i<gdata->n_ele;i++) {
00122 ierr = ProgressBar(i,gdata->n_ele,10);
00123
00124 int t_nnb;
00125 t_nnb=ia[i+1]-ia[i];
00126
00127
00128 for (int j=0; j<NV; j++) {
00129 t_elevs[j]=gdata->elevert[NV*i+j];
00130
00131
00132
00133
00134
00135 }
00136
00137
00138 ierr = PetscSortInt(NV,(PetscInt*)t_elevs);CHKERRQ(ierr);
00139
00140
00141 t_facvert[0] = t_elevs[0];
00142 t_facvert[1] = t_elevs[1];
00143 t_facvert[2] = t_elevs[2];
00144 t_facvert[3] = t_elevs[3];
00145
00146 t_facvert[4] = t_elevs[0];
00147 t_facvert[5] = t_elevs[1];
00148 t_facvert[6] = t_elevs[3];
00149 t_facvert[7] = t_elevs[2];
00150
00151 t_facvert[8] = t_elevs[0];
00152 t_facvert[9] = t_elevs[2];
00153 t_facvert[10] = t_elevs[3];
00154 t_facvert[11] = t_elevs[1];
00155
00156 t_facvert[12] = t_elevs[1];
00157 t_facvert[13] = t_elevs[2];
00158 t_facvert[14] = t_elevs[3];
00159 t_facvert[15] = t_elevs[0];
00160
00161
00162 for (int j=0; j<NF; j++) {
00163 int t_found;
00164 t_found=0;
00165
00166
00167 for (int k=0; k<t_nnb; k++) {
00168
00169 int t_nbid;
00170 t_nbid=ja[ia[i]+k];
00171
00172
00173 int x,vshared[NV];
00174 x=tetnb(gdata->elevert+NV*i,gdata->elevert+NV*t_nbid,vshared);
00175
00176 assert(x==NN);
00177
00178 if ((t_facvert[NV*j+0]==vshared[0]) &&
00179 (t_facvert[NV*j+1]==vshared[1]) &&
00180 (t_facvert[NV*j+2]==vshared[2])
00181 ) {
00182
00183 t_found=1;
00184
00185
00186
00187 if (t_nbid<i) {
00188
00189 for (int m=0; m<NF; m++) {
00190 if (facnb[elefac[t_nbid*NF+m]*2+1]==i) {
00191 elefac[i*NF+j]=elefac[t_nbid*NF+m];
00192 break;
00193 }
00194 }
00195 }
00196
00197 else {
00198 facvert[t_facmax*NN+0]=vshared[0];
00199 facvert[t_facmax*NN+1]=vshared[1];
00200 facvert[t_facmax*NN+2]=vshared[2];
00201 facnb[t_facmax*2+0]=i;
00202 facnb[t_facmax*2+1]=t_nbid;
00203 elefac[i*NF+j]=t_facmax;
00204 t_facmax++;
00205 }
00206 break;
00207 }
00208 if (t_found) break;
00209 }
00210
00211 if (!t_found) {
00212
00213
00214
00215
00216
00217
00218
00219 assert(t_nnb<NF);
00220
00221 PetscReal *x1,*x2,*x3,*x4;
00222 x1=gdata->vertxyz+(t_facvert[NV*j+0])*ND;
00223 x2=gdata->vertxyz+(t_facvert[NV*j+1])*ND;
00224 x3=gdata->vertxyz+(t_facvert[NV*j+2])*ND;
00225 x4=gdata->vertxyz+(t_facvert[NV*j+3])*ND;
00226
00227 facvert[t_facmax*NN+0]=t_facvert[NV*j+0];
00228
00229 if (tetvol(x1,x2,x3,x4)<0.0) {
00230 facvert[t_facmax*NN+1]=t_facvert[NV*j+1];
00231 facvert[t_facmax*NN+2]=t_facvert[NV*j+2];
00232 }
00233
00234 else {
00235 facvert[t_facmax*NN+1]=t_facvert[NV*j+2];
00236 facvert[t_facmax*NN+2]=t_facvert[NV*j+1];
00237 }
00238 facnb[t_facmax*2+0]=i;
00239 facnb[t_facmax*2+1]=C_BND;
00240 elefac[i*NF+j]=t_facmax;
00241 t_facmax++;
00242 }
00243 }
00244 }
00245 ierr = ProgressBar(1,1,10);
00246
00247 PetscLogDouble t_t3;
00248 ierr = PetscGetTime(&t_t3);CHKERRQ(ierr);
00249 PetscPrintf(PETSC_COMM_WORLD,
00250 "Search for surface faces and nodes took %g s = %g min\n",
00251 t_t3-t_t2,
00252 (t_t3-t_t2)/60.0
00253 );
00254
00255 #ifdef ADDONS
00256 PetscTruth flg;
00257 int mergepoints;
00258 ierr = PetscOptionsGetInt(PETSC_NULL,"-hdemag_mergepoints",(PetscInt*)&mergepoints,&flg);CHKERRQ(ierr);
00259 if (flg!=PETSC_TRUE) {
00260 mergepoints=1;
00261 PetscPrintf(PETSC_COMM_WORLD,
00262 "Option -hdemag_mergepoints not found, using default value: %i\n",
00263 mergepoints
00264 );
00265 }
00266 if (gdata->vertcoupl==PETSC_NULL) {
00267 mergepoints=0;
00268 PetscPrintf(PETSC_COMM_WORLD,"No coupled vertices found, setting -hdemag_mergepoints %i!\n",mergepoints);
00269 }
00270
00271 int t_cnt;
00272 t_cnt=0;
00273 if (mergepoints) {
00274
00275 for (int i=0;i<n_fac;i++) {
00276 for (int j=0;j<NN;j++) {
00277 if (gdata->vertcoupl[facvert[NN*i+j]]>=0) {
00278
00279 facvert[NN*i+j]=PetscMin(facvert[NN*i+j],gdata->vertcoupl[facvert[NN*i+j]]);
00280 }
00281 }
00282 #define swap(a,b) {int tswap;tswap=a;a=b;b=tswap;}
00283 if (facvert[NN*i+0]>facvert[NN*i+1]) {
00284 swap(facvert[NN*i+0],facvert[NN*i+1]);
00285 swap(facvert[NN*i+1],facvert[NN*i+2]);
00286 }
00287 if (facvert[NN*i+0]>facvert[NN*i+2]) {
00288 swap(facvert[NN*i+0],facvert[NN*i+1]);
00289 swap(facvert[NN*i+1],facvert[NN*i+2]);
00290 }
00291 if (facvert[NN*i+0]>facvert[NN*i+1]) {
00292 swap(facvert[NN*i+0],facvert[NN*i+1]);
00293 swap(facvert[NN*i+1],facvert[NN*i+2]);
00294 }
00295 }
00296
00297 PetscPrintf(PETSC_COMM_WORLD,"Searching for coupled faces/nodes...\n");
00298
00299 int nsfac=0;
00300 int *sfac;
00301 ierr = PetscMalloc(n_fac*sizeof(int),&sfac);CHKERRQ(ierr);
00302
00303 for (int i=0;i<n_fac;i++) {
00304 if (facnb[i*2+1]!=C_BND ||
00305 gdata->vertcoupl[facvert[NN*i+0]]==-1 ||
00306 gdata->vertcoupl[facvert[NN*i+1]]==-1 ||
00307 gdata->vertcoupl[facvert[NN*i+2]]==-1) continue;
00308
00309 sfac[nsfac]=i;
00310 nsfac++;
00311 }
00312
00313 #ifdef USEQSORT
00314 PetscPrintf(PETSC_COMM_WORLD,"Sorting surface triangles...\n");
00315 qsort(sfac,nsfac,sizeof(int),CompareFaces);
00316 #endif
00317
00318 PetscPrintf(PETSC_COMM_WORLD,"Checking %i coupled faces/nodes...\n",nsfac);
00319
00320 for (int i2=0;i2<nsfac-1;i2++) {
00321 ierr = ProgressBar(i2,nsfac-1,10);
00322
00323 int i;
00324 i=sfac[i2];
00325
00326
00327 for (int k2=i2+1;k2<nsfac;k2++) {
00328
00329 int k;
00330 k=sfac[k2];
00331
00332 #ifdef USEQSORT
00333 if (facvert[NN*k+0]>facvert[NN*i+0]) break;
00334 #endif
00335
00336
00337
00338
00339
00340
00341 if (
00342 facvert[NN*i+0]==facvert[NN*k+0] &&
00343 facvert[NN*i+1]==facvert[NN*k+2] &&
00344 facvert[NN*i+2]==facvert[NN*k+1]
00345 ) {
00346
00347 facnb[i*2+1]=facnb[k*2+0];
00348 facnb[k*2+1]=facnb[i*2+0];
00349
00350 t_cnt+=2;
00351 break;
00352 }
00353 }
00354 }
00355 ierr = ProgressBar(1,1,10);
00356
00357 ierr = PetscFree(sfac);CHKERRQ(ierr);
00358
00359 PetscPrintf(PETSC_COMM_WORLD,"Retagged %i coupled surfaces\n",t_cnt);
00360 }
00361
00362 t_t2=t_t3;
00363 ierr = PetscGetTime(&t_t3);CHKERRQ(ierr);
00364 PetscPrintf(PETSC_COMM_WORLD,
00365 "Search for coupled faces/nodes took %g s = %g min\n",
00366 t_t3-t_t2,
00367 (t_t3-t_t2)/60.0
00368 );
00369 #endif
00370
00371 ierr = PetscMalloc(gdata->n_vert*sizeof(int),&gdata->vertbndg2bnd);CHKERRQ(ierr);
00372
00373
00374 for (int i=0;i<gdata->n_vert;i++) {
00375 gdata->vertbndg2bnd[i]=C_INT;
00376 }
00377
00378
00379 gdata->n_vert_bnd=0;
00380 for (int i=0;i<n_fac;i++) {
00381 if (facnb[i*2+1]==C_BND) {
00382 for (int j=0;j<NN;j++) {
00383 int k;
00384 k=facvert[NN*i+j];
00385 #ifdef ADDONS
00386 if (mergepoints && gdata->vertcoupl[k]>=0) {
00387
00388 assert(k<gdata->vertcoupl[k]);
00389 }
00390 #endif
00391
00392 if (gdata->vertbndg2bnd[k]==C_INT) {
00393 gdata->vertbndg2bnd[k]=C_BND;
00394 gdata->n_vert_bnd++;
00395 }
00396 }
00397 }
00398 }
00399 PetscPrintf(PETSC_COMM_WORLD,
00400 "Boundary matrix: %i vertices, %g matrix elements, %i MB \n",
00401 gdata->n_vert_bnd,
00402 1.0*gdata->n_vert_bnd*gdata->n_vert_bnd,
00403 gdata->n_vert_bnd/1024*gdata->n_vert_bnd/1024*sizeof(PetscReal)+1
00404 );
00405
00406
00407 gdata->n_bnd_fac=0;
00408 for (int i=0;i<gdata->n_ele;i++) {
00409 for (int j=0;j<NF;j++) {
00410 if (facnb[elefac[i*NF+j]*2+1]==C_BND) {
00411 gdata->n_bnd_fac++;
00412 }
00413 }
00414 }
00415 gdata->ln_bnd_fac=gdata->n_bnd_fac;
00416 PetscPrintf(PETSC_COMM_WORLD,"Number of boundary faces: %i\n",gdata->n_bnd_fac);
00417
00418
00419 ierr = PetscMalloc(gdata->n_bnd_fac*NN*sizeof(int),&gdata->bndfacvert);CHKERRQ(ierr);
00420 int t_cnt2;
00421 t_cnt2=0;
00422 for (int i=0;i<n_fac;i++) {
00423 if (facnb[i*2+1]==C_BND) {
00424 for (int j=0;j<NN;j++) {
00425 gdata->bndfacvert[t_cnt2*NN+j]=facvert[i*NN+j];
00426 }
00427 t_cnt2++;
00428 }
00429 }
00430 assert(t_cnt2==gdata->n_bnd_fac);
00431
00432
00433 int t_cnt3;
00434 t_cnt3=0;
00435 for (int i=0; i<gdata->n_vert; i++) {
00436 if (gdata->vertbndg2bnd[i]==C_BND) {
00437 gdata->vertbndg2bnd[i]=t_cnt3++;
00438 }
00439 }
00440 assert(t_cnt3==gdata->n_vert_bnd);
00441
00442 ierr = PetscFree(elefac);CHKERRQ(ierr);
00443 ierr = PetscFree(facnb);CHKERRQ(ierr);
00444 ierr = PetscFree(facvert);CHKERRQ(ierr);
00445
00446 ierr = PetscFree(ia);CHKERRQ(ierr);
00447 ierr = PetscFree(ja);CHKERRQ(ierr);
00448
00449
00450 MagparFunctionLogReturn(0);
00451 }