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: elevertvol.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 PetscReal tetgrad(PetscReal *x1,PetscReal *x2,PetscReal *x3,PetscReal *x4,PetscReal D_etaj[][ND])
00031 {
00032 PetscReal t_elevol;
00033
00034 t_elevol=tetvol(x1,x2,x3,x4);
00035
00036
00037 D_etaj[0][0]=-(x3[1]*x2[2]-x4[1]*x2[2]-x2[1]*x3[2]+x4[1]*x3[2]+x2[1]*x4[2]-x3[1]*x4[2]);
00038 D_etaj[0][1]=+(x3[0]*x2[2]-x4[0]*x2[2]-x2[0]*x3[2]+x4[0]*x3[2]+x2[0]*x4[2]-x3[0]*x4[2]);
00039 D_etaj[0][2]=-(x3[0]*x2[1]-x4[0]*x2[1]-x2[0]*x3[1]+x4[0]*x3[1]+x2[0]*x4[1]-x3[0]*x4[1]);
00040
00041 D_etaj[1][0]=+(x4[1]*x3[2]-x1[1]*x3[2]-x3[1]*x4[2]+x1[1]*x4[2]+x3[1]*x1[2]-x4[1]*x1[2]);
00042 D_etaj[1][1]=-(x4[0]*x3[2]-x1[0]*x3[2]-x3[0]*x4[2]+x1[0]*x4[2]+x3[0]*x1[2]-x4[0]*x1[2]);
00043 D_etaj[1][2]=+(x4[0]*x3[1]-x1[0]*x3[1]-x3[0]*x4[1]+x1[0]*x4[1]+x3[0]*x1[1]-x4[0]*x1[1]);
00044
00045 D_etaj[2][0]=-(x1[1]*x4[2]-x2[1]*x4[2]-x4[1]*x1[2]+x2[1]*x1[2]+x4[1]*x2[2]-x1[1]*x2[2]);
00046 D_etaj[2][1]=+(x1[0]*x4[2]-x2[0]*x4[2]-x4[0]*x1[2]+x2[0]*x1[2]+x4[0]*x2[2]-x1[0]*x2[2]);
00047 D_etaj[2][2]=-(x1[0]*x4[1]-x2[0]*x4[1]-x4[0]*x1[1]+x2[0]*x1[1]+x4[0]*x2[1]-x1[0]*x2[1]);
00048
00049 D_etaj[3][0]=+(x2[1]*x1[2]-x3[1]*x1[2]-x1[1]*x2[2]+x3[1]*x2[2]+x1[1]*x3[2]-x2[1]*x3[2]);
00050 D_etaj[3][1]=-(x2[0]*x1[2]-x3[0]*x1[2]-x1[0]*x2[2]+x3[0]*x2[2]+x1[0]*x3[2]-x2[0]*x3[2]);
00051 D_etaj[3][2]=+(x2[0]*x1[1]-x3[0]*x1[1]-x1[0]*x2[1]+x3[0]*x2[1]+x1[0]*x3[1]-x2[0]*x3[1]);
00052
00053 for (int i=0;i<NV;i++) {
00054 for (int j=0;j<ND;j++) {
00055 D_etaj[i][j] /= 6.0*t_elevol;
00056 }
00057 }
00058
00059
00060 for (int j=0;j<ND;j++) {
00061 if (PetscAbsReal(D_etaj[0][j]+D_etaj[1][j]+D_etaj[2][j]+D_etaj[3][j])>=D_EPS*1000) {
00062 PetscPrintf(PETSC_COMM_WORLD,
00063 "etaj: %g %g %g %g sum= %g\n",
00064 D_etaj[0][j],D_etaj[1][j],D_etaj[2][j],D_etaj[3][j],
00065 D_etaj[0][j]+D_etaj[1][j]+D_etaj[2][j]+D_etaj[3][j]
00066 );
00067 }
00068 }
00069
00070 return(0);
00071 }
00072
00073
00074 PetscReal tetqual(PetscReal *x1,PetscReal *x2,PetscReal *x3,PetscReal *x4)
00075 {
00076 PetscReal result;
00077
00078
00079
00080
00081
00082
00083
00084
00085
00086 result=tetvol(x1,x2,x3,x4)/(sqrt(2.0)/12.0*
00087 pow(
00088 (
00089 sqrt((x1[0]-x2[0])*(x1[0]-x2[0])+(x1[1]-x2[1])*(x1[1]-x2[1])+(x1[2]-x2[2])*(x1[2]-x2[2]))+
00090 sqrt((x1[0]-x3[0])*(x1[0]-x3[0])+(x1[1]-x3[1])*(x1[1]-x3[1])+(x1[2]-x3[2])*(x1[2]-x3[2]))+
00091 sqrt((x1[0]-x4[0])*(x1[0]-x4[0])+(x1[1]-x4[1])*(x1[1]-x4[1])+(x1[2]-x4[2])*(x1[2]-x4[2]))+
00092 sqrt((x2[0]-x3[0])*(x2[0]-x3[0])+(x2[1]-x3[1])*(x2[1]-x3[1])+(x2[2]-x3[2])*(x2[2]-x3[2]))+
00093 sqrt((x2[0]-x4[0])*(x2[0]-x4[0])+(x2[1]-x4[1])*(x2[1]-x4[1])+(x2[2]-x4[2])*(x2[2]-x4[2]))+
00094 sqrt((x3[0]-x4[0])*(x3[0]-x4[0])+(x3[1]-x4[1])*(x3[1]-x4[1])+(x3[2]-x4[2])*(x3[2]-x4[2]))
00095 )/6.0,
00096 3.0
00097 )
00098 );
00099
00100 return(result);
00101 }
00102
00103
00104 PetscReal edgestat(int nele,int *ele,PetscReal *xyz,PetscReal *min,PetscReal *max,PetscReal *avg)
00105 {
00106 PetscReal *x1,*x2,*x3,*x4;
00107 PetscReal edge[ND];
00108 PetscReal elen;
00109 PetscReal elenavg=0;
00110 PetscReal elenmin=PETSC_MAX;
00111 PetscReal elenmax=PETSC_MIN;
00112
00113 MagparFunctionInfoBegin;
00114
00115 for (int i=0; i<nele; i++) {
00116
00117
00118 int *t_v;
00119 t_v=ele+i*NV;
00120
00121 x1=xyz+t_v[0]*ND;
00122 x2=xyz+t_v[1]*ND;
00123 x3=xyz+t_v[2]*ND;
00124 x4=xyz+t_v[3]*ND;
00125
00126
00127 my_dcopy(ND,x1,1,edge,1);
00128 my_daxpy(ND,-1.0,x2,1,edge,1);
00129 elen=my_dnrm2(ND,edge,1);
00130 elenmax=PetscMax(elenmax,elen);
00131 elenmin=PetscMin(elenmin,elen);
00132 elenavg += elen;
00133
00134 my_dcopy(ND,x1,1,edge,1);
00135 my_daxpy(ND,-1.0,x3,1,edge,1);
00136 elen=my_dnrm2(ND,edge,1);
00137 elenmax=PetscMax(elenmax,elen);
00138 elenmin=PetscMin(elenmin,elen);
00139 elenavg += elen;
00140
00141 my_dcopy(ND,x1,1,edge,1);
00142 my_daxpy(ND,-1.0,x4,1,edge,1);
00143 elen=my_dnrm2(ND,edge,1);
00144 elenmax=PetscMax(elenmax,elen);
00145 elenmin=PetscMin(elenmin,elen);
00146 elenavg += elen;
00147
00148 my_dcopy(ND,x2,1,edge,1);
00149 my_daxpy(ND,-1.0,x3,1,edge,1);
00150 elen=my_dnrm2(ND,edge,1);
00151 elenmax=PetscMax(elenmax,elen);
00152 elenmin=PetscMin(elenmin,elen);
00153 elenavg += elen;
00154
00155 my_dcopy(ND,x2,1,edge,1);
00156 my_daxpy(ND,-1.0,x4,1,edge,1);
00157 elen=my_dnrm2(ND,edge,1);
00158 elenmax=PetscMax(elenmax,elen);
00159 elenmin=PetscMin(elenmin,elen);
00160 elenavg += elen;
00161
00162 my_dcopy(ND,x3,1,edge,1);
00163 my_daxpy(ND,-1.0,x4,1,edge,1);
00164 elen=my_dnrm2(ND,edge,1);
00165 elenmax=PetscMax(elenmax,elen);
00166 elenmin=PetscMin(elenmin,elen);
00167 elenavg += elen;
00168 }
00169
00170 ierr = PetscGlobalMin(&elenmin,min,PETSC_COMM_WORLD);CHKERRQ(ierr);
00171 ierr = PetscGlobalMax(&elenmax,max,PETSC_COMM_WORLD);CHKERRQ(ierr);
00172 ierr = PetscGlobalSum(&elenavg,avg,PETSC_COMM_WORLD);CHKERRQ(ierr);
00173
00174 MagparFunctionInfoReturn(0);
00175 }
00176
00177
00178 int EleVertVol(GridData *gdata)
00179 {
00180 MagparFunctionLogBegin;
00181
00182 int rank,size;
00183 ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRQ(ierr);
00184 ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRQ(ierr);
00185
00186 PetscReal *vprop;
00187 ierr = PetscMalloc(gdata->n_prop*sizeof(PetscReal),&vprop);CHKERRQ(ierr);
00188 for (int i=0;i<gdata->n_prop;i++) {
00189 vprop[i]=0.0;
00190 }
00191
00192 ierr = VecCreate(PETSC_COMM_WORLD,&gdata->elevol);CHKERRQ(ierr);
00193 ierr = VecSetSizes(gdata->elevol,gdata->ln_ele,gdata->n_ele);CHKERRQ(ierr);
00194 ierr = VecSetFromOptions(gdata->elevol);CHKERRQ(ierr);
00195
00196 Vec elequal;
00197 ierr = VecDuplicate(gdata->elevol,&elequal);CHKERRQ(ierr);
00198
00199 ierr = VecCreate(PETSC_COMM_WORLD,&gdata->vertvol);CHKERRQ(ierr);
00200 ierr = VecSetSizes(gdata->vertvol,gdata->ln_vert,gdata->n_vert);CHKERRQ(ierr);
00201 ierr = VecSetFromOptions(gdata->vertvol);CHKERRQ(ierr);
00202
00203 ierr = VecDuplicate(gdata->M,&gdata->VMs3);CHKERRQ(ierr);
00204
00205 PetscReal Javg=0.0;
00206 PetscReal Vmag=0.0;
00207 for (int i=0; i<gdata->ln_ele; i++) {
00208
00209 int *t_v;
00210 t_v=gdata->elevert+i*NV;
00211
00212 PetscReal *x1,*x2,*x3,*x4;
00213 x1=gdata->vertxyz+t_v[0]*ND;
00214 x2=gdata->vertxyz+t_v[1]*ND;
00215 x3=gdata->vertxyz+t_v[2]*ND;
00216 x4=gdata->vertxyz+t_v[3]*ND;
00217
00218
00219 PetscReal T;
00220 T=tetvol(x1,x2,x3,x4);
00221
00222
00223 if (T<0.0) {
00224
00225
00226
00227
00228
00229
00230 int j;
00231 j=t_v[0];
00232 t_v[0]=t_v[1];
00233 t_v[1]=j;
00234 T=-T;
00235
00236 x1=gdata->vertxyz+t_v[0]*ND;
00237 x2=gdata->vertxyz+t_v[1]*ND;
00238 }
00239
00240 if (gdata->propdat[NP*gdata->eleprop[i]+4]>0.0) {
00241 Vmag += T;
00242 }
00243
00244 ierr = VecSetValue(
00245 gdata->elevol,
00246 gdata->elel2g[i],
00247 T,
00248 INSERT_VALUES
00249 );CHKERRQ(ierr);
00250
00251 for (int j=0;j<NV;j++) {
00252 ierr = VecSetValue(
00253 gdata->vertvol,
00254 t_v[j],
00255 T/NV,
00256 ADD_VALUES
00257 );CHKERRQ(ierr);
00258
00259 PetscReal matele;
00260 matele=T/NV*gdata->propdat[NP*gdata->eleprop[i]+4];
00261 if (PetscAbsReal(matele)>D_EPS) {
00262 for (int k=0;k<ND;k++) {
00263 ierr = VecSetValue(
00264 gdata->VMs3,
00265 gdata->elevert[NV*i+j]*ND+k,
00266 matele,
00267 ADD_VALUES
00268 );CHKERRQ(ierr);
00269 }
00270 }
00271 }
00272 Javg += T*gdata->propdat[NP*gdata->eleprop[i]+4]*gdata->hscale;
00273 vprop[gdata->eleprop[i]] += T;
00274
00275 PetscReal t_qual;
00276 t_qual=tetqual(x1,x2,x3,x4);
00277 ierr = VecSetValue(
00278 elequal,
00279 gdata->elel2g[i],
00280 t_qual,
00281 INSERT_VALUES
00282 );CHKERRQ(ierr);
00283 }
00284
00285 ierr = VecAssemblyBegin(gdata->elevol);CHKERRQ(ierr);
00286 ierr = VecAssemblyEnd(gdata->elevol);CHKERRQ(ierr);
00287 ierr = VecAssemblyBegin(gdata->vertvol);CHKERRQ(ierr);
00288 ierr = VecAssemblyEnd(gdata->vertvol);CHKERRQ(ierr);
00289 ierr = VecAssemblyBegin(gdata->VMs3);CHKERRQ(ierr);
00290 ierr = VecAssemblyEnd(gdata->VMs3);CHKERRQ(ierr);
00291
00292 for (int i=0;i<gdata->n_prop;i++) {
00293 PetscReal temp;
00294 ierr = PetscGlobalSum(&vprop[i],&temp,PETSC_COMM_WORLD);CHKERRQ(ierr);
00295 vprop[i]=temp;
00296 }
00297
00298
00299
00300 char fmesh[256];
00301 FILE *fd;
00302 ierr = PetscSNPrintf(fmesh,255,"%s.%04i.%s",gdata->simname,gdata->inp,"felog");CHKERRQ(ierr);
00303 ierr = PetscFOpen(PETSC_COMM_WORLD,fmesh,"w",&fd);CHKERRQ(ierr);
00304 if (!rank && !fd) {
00305 SETERRQ1(PETSC_ERR_FILE_OPEN,"Cannot open %s\n",fmesh);
00306 }
00307 PetscPrintf(PETSC_COMM_WORLD,"Opening file %s\n",fmesh);
00308
00309
00310
00311 ierr = PetscFPrintf(PETSC_COMM_WORLD,fd,
00312 "simname: %s\n"
00313 "n_vert: %7d\n"
00314 "n_ele: %7d\n"
00315 "n_vert_bnd: %7d\n"
00316 "n_bnd_fac: %7d\n"
00317 "n_prop: %7d\n",
00318 gdata->simname,
00319 gdata->n_vert,
00320 gdata->n_ele,
00321 gdata->n_vert_bnd,
00322 gdata->n_bnd_fac,
00323 gdata->n_prop
00324 );
00325
00326
00327 PetscReal bbox[2*ND];
00328 int pid[1]={-1};
00329
00330 ierr = calcbbox(gdata->ln_ele,gdata->elevert,PETSC_NULL,gdata->vertxyz,1,pid,bbox);CHKERRQ(ierr);
00331
00332 ierr = PetscFPrintf(PETSC_COMM_WORLD,fd,
00333 "\nbounding box: (%g,%g,%g) (%g,%g,%g)\n",
00334 bbox[0],bbox[1],bbox[2],
00335 bbox[ND+0],bbox[ND+1],bbox[ND+2]
00336 );CHKERRQ(ierr);
00337
00338
00339
00340 PetscReal Vtot,Vmin,Vmax;
00341 int t_idxmin,t_idxmax;
00342 ierr = VecMin(gdata->elevol,(PetscInt*)&t_idxmin,&Vmin);CHKERRQ(ierr);
00343 ierr = VecMax(gdata->elevol,(PetscInt*)&t_idxmax,&Vmax);CHKERRQ(ierr);
00344 ierr = VecSum(gdata->elevol,&Vtot);CHKERRQ(ierr);
00345 gdata->totvol=Vtot;
00346
00347 ierr = PetscFPrintf(PETSC_COMM_WORLD,fd,
00348 "\nelevol:\n"
00349 " id_min: %i\n"
00350 " id_max: %i\n"
00351 " Vmax: %g\n"
00352 " Vmin: %g\n"
00353 " ratio: %g\n"
00354 " Vavg: %g\n"
00355 " Vtot: %g\n",
00356 t_idxmin,t_idxmax,Vmax,Vmin,Vmax/Vmin,Vtot/gdata->n_ele,Vtot
00357 );
00358
00359
00360
00361 ierr = VecMin(gdata->vertvol,(PetscInt*)&t_idxmin,&Vmin);CHKERRQ(ierr);
00362 ierr = VecMax(gdata->vertvol,(PetscInt*)&t_idxmax,&Vmax);CHKERRQ(ierr);
00363 ierr = VecSum(gdata->vertvol,&Vtot);CHKERRQ(ierr);
00364
00365
00366 assert(PetscAbsReal(gdata->totvol-Vtot) < Vtot*1e-7);
00367
00368 ierr = PetscFPrintf(PETSC_COMM_WORLD,fd,
00369 "\nvertvol:\n"
00370 " id_min: %i\n"
00371 " id_max: %i\n"
00372 " Vmax: %g\n"
00373 " Vmin: %g\n"
00374 " ratio: %g\n"
00375 " Vavg: %g\n"
00376 " Vtot: %g\n",
00377 t_idxmin,t_idxmax,Vmax,Vmin,Vmax/Vmin,Vtot/gdata->n_ele,Vtot
00378 );
00379
00380
00381
00382 PetscReal elenmax,elenmin,elenavg;
00383 edgestat(gdata->ln_ele,gdata->elevert,gdata->vertxyz,&elenmin,&elenmax,&elenavg);
00384 elenavg /= 6*gdata->n_ele;
00385 ierr = PetscFPrintf(PETSC_COMM_WORLD,fd,
00386 "\n"
00387 "edge_len_min: %g\n"
00388 "edge_len_max: %g\n"
00389 "edge_len_avg: %g\n",
00390 elenmin,elenmax,elenavg
00391 );
00392 gdata->elenmax=elenmax;
00393
00394
00395
00396
00397 ierr = VecMin(elequal,(PetscInt*)&t_idxmin,&Vmin);CHKERRQ(ierr);
00398 ierr = VecMax(elequal,(PetscInt*)&t_idxmax,&Vmax);CHKERRQ(ierr);
00399 ierr = VecSum(elequal,&Vtot);CHKERRQ(ierr);
00400 ierr = PetscFPrintf(PETSC_COMM_WORLD,fd,
00401 "\nelequal:\n"
00402 " min: %g\n"
00403 " max: %g\n"
00404 " avg: %g\n",
00405 Vmin,Vmax,Vtot/gdata->n_ele
00406 );
00407
00408 ierr = VecDestroy(elequal);CHKERRQ(ierr);
00409
00410
00411
00412 Vtot=0.0;
00413 ierr = PetscFPrintf(PETSC_COMM_WORLD,fd,
00414 "\nvolume by property id\n\n"
00415 " pid | vol\n"
00416 "-----------\n"
00417 );
00418 for (int i=0;i<gdata->n_prop;i++) {
00419 ierr = PetscFPrintf(PETSC_COMM_WORLD,fd,
00420 "%5d %g\n",
00421 i+1,vprop[i]
00422 );
00423 Vtot += vprop[i];
00424 }
00425 ierr = PetscFPrintf(PETSC_COMM_WORLD,fd,
00426 "-----------\n"
00427 " sum = %g\n",
00428 Vtot
00429 );
00430 #ifdef PETSC_USE_SINGLE
00431 assert(PetscAbsReal(gdata->totvol-Vtot) < Vtot*1e-3);
00432 #else
00433 assert(PetscAbsReal(gdata->totvol-Vtot) < Vtot*1e-7);
00434 #endif
00435
00436
00437
00438 ierr = PetscGlobalSum(&Vmag,&Vtot,PETSC_COMM_WORLD);CHKERRQ(ierr);
00439 gdata->totvol=Vtot;
00440
00441
00442
00443 PetscReal sum;
00444 ierr = PetscGlobalSum(&Javg,&sum,PETSC_COMM_WORLD);CHKERRQ(ierr);
00445 Javg=sum/Vtot;
00446 ierr = PetscFPrintf(PETSC_COMM_WORLD,fd,
00447 "\naverage magnetization: %g T\n",
00448 Javg
00449 );
00450
00451
00452
00453 ierr = PetscFClose(PETSC_COMM_WORLD,fd);CHKERRQ(ierr);
00454 ierr = PetscFree(vprop);CHKERRQ(ierr);
00455
00456 MagparFunctionLogReturn(0);
00457 }
00458