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: writedatadat.c 3051 2010-04-14 17:58:54Z scholz $\n\n";
00025 static char const Td[] = "$Today: " __FILE__ " " __DATE__ " " __TIME__ " $\n\n";
00026
00027 #include "magpario.h"
00028 #include "util/util.h"
00029 #ifdef ADDONS
00030 #include "addons/addons.h"
00031 #endif
00032
00033 #define NVDATA 18
00034 #define NEDATA 1
00035 #define NEDATA2 0
00036
00037 static int doinit=1;
00038 static Mat AIpolline;
00039 static Vec VIpolline;
00042 int WriteDatInit(GridData *gdata)
00043 {
00044 int i,j,k,l;
00045
00046 int nmax;
00047 PetscTruth flg;
00048
00049 int res;
00050 PetscReal *vertxyz2;
00051 PetscReal line_p[ND], line_v[ND];
00052
00053 PetscReal axes2[ND*ND];
00054
00055 int cnt,*sliceele;
00056 PetscReal sgnx,sgny;
00057 PetscReal x2min[ND]={0,0,PETSC_MAX};
00058 PetscReal x2max[ND]={0,0,PETSC_MIN};
00059 PetscReal xmin[ND],xmax[ND];
00060 PetscReal line_length;
00061 PetscReal *d,p[ND];
00062 PetscReal rhs[ND];
00063 int t_ele;
00064
00065 PetscReal matele;
00066
00067
00068 MagparFunctionLogBegin;
00069
00070 int rank,size;
00071 ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRQ(ierr);
00072 ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRQ(ierr);
00073
00074 ierr = PetscMalloc(ND*gdata->n_vert*sizeof(PetscReal),&vertxyz2);CHKERRQ(ierr);
00075
00076 nmax=ND;
00077 ierr = PetscOptionsGetRealArray(PETSC_NULL,"-line_p",line_p,(PetscInt*)&nmax,&flg);CHKERRQ(ierr);
00078 if (flg!=PETSC_TRUE) {
00079 line_p[0]=line_p[1]=line_p[2]=PETSC_MAX;
00080 PetscPrintf(PETSC_COMM_WORLD,
00081 "Option -line_p not found, using default value: (%g,%g,%g)\n",
00082 line_p[0],line_p[1],line_p[2]
00083 );
00084 nmax=ND;
00085 }
00086 if (nmax!=ND)
00087 SETERRQ(PETSC_ERR_ARG_INCOMP,"Wrong number of values for option -line_p!\n");
00088
00089 nmax=ND;
00090 ierr = PetscOptionsGetRealArray(PETSC_NULL,"-line_v",line_v,(PetscInt*)&nmax,&flg);CHKERRQ(ierr);
00091 if (flg!=PETSC_TRUE) {
00092 line_v[0]=1;line_v[1]=line_v[2]=0;
00093 PetscPrintf(PETSC_COMM_WORLD,
00094 "Option -line_v not found, using default value: (%g,%g,%g)\n",
00095 line_v[0],line_v[1],line_v[2]
00096 );
00097 nmax=ND;
00098 }
00099 if (nmax!=ND)
00100 SETERRQ(PETSC_ERR_ARG_INCOMP,"Wrong number of values for option -line_v!\n");
00101
00102 PetscPrintf(PETSC_COMM_WORLD,
00103 "sampling line: (x,y,z)=(%g,%g,%g)+n*(%g,%g,%g)\n",
00104 line_p[0],line_p[1],line_p[2],
00105 line_v[0],line_v[1],line_v[2]
00106 );
00107
00108
00109 PetscReal nrm;
00110 nrm=1.0/my_dnrm2(ND,line_v,1);
00111 my_dscal(ND,nrm,line_v,1);
00112
00113
00114 PetscPrintf(PETSC_COMM_WORLD,
00115 "normalized sampling line vector: (%g,%g,%g)\n",
00116 line_v[0],line_v[1],line_v[2]
00117 );
00118
00119 PetscPrintf(PETSC_COMM_WORLD,"shift and tilt to avoid cutting 'between' elements\n");
00120 #define DIST 1e-12
00121 line_p[0]+=DIST;line_v[0]+=DIST;
00122 line_p[1]+=DIST;line_v[1]+=DIST;
00123 line_p[2]+=DIST;line_v[2]+=DIST;
00124 my_dscal(ND,1.0+DIST,line_v,1);
00125 my_dscal(ND,1.0+DIST,line_p,1);
00126
00127
00128 AxesRotation(line_p,line_v,axes2,0,gdata->n_vert,gdata->vertxyz,vertxyz2);
00129
00130 ierr = PetscMalloc(gdata->ln_ele*sizeof(int),&sliceele);CHKERRQ(ierr);
00131
00132 cnt=0;
00133 for (i=0;i<gdata->ln_ele;i++) {
00134
00135 sgnx=vertxyz2[ND*gdata->elevert[NV*i+0]+0];
00136 sgny=vertxyz2[ND*gdata->elevert[NV*i+0]+1];
00137
00138 for (j=1;j<NV;j++) {
00139 if (sgnx*vertxyz2[ND*gdata->elevert[NV*i+j]+0] <= 0.0) break;
00140 }
00141 if (j<NV) {
00142 for (j=1;j<NV;j++) {
00143
00144
00145
00146
00147 if (sgny*vertxyz2[ND*gdata->elevert[NV*i+j]+1] <= 0.0) break;
00148 }
00149
00150 if (j<NV) {
00151
00152 sliceele[cnt++]=i;
00153 }
00154 }
00155 }
00156 PetscSynchronizedPrintf(PETSC_COMM_WORLD,
00157 "<%i>Found %i elements cut by sampling line\n",
00158 rank,cnt
00159 );
00160 PetscSynchronizedFlush(PETSC_COMM_WORLD);
00161
00162 int sum;
00163 ierr=MPI_Allreduce((int*)&cnt,(int*)&sum,1,MPI_INT,MPI_SUM,PETSC_COMM_WORLD);CHKERRQ(ierr);
00164 if (sum==0) {
00165 AIpolline=PETSC_NULL;
00166
00167 PetscPrintf(PETSC_COMM_WORLD,
00168 "Sampling line outside the FE mesh.\nNo *.d files will be stored!\n"
00169 );
00170
00171 ierr = PetscFree(sliceele);CHKERRQ(ierr);
00172 ierr = PetscFree(vertxyz2);CHKERRQ(ierr);
00173
00174 MagparFunctionLogReturn(0);
00175 }
00176
00177
00178 ierr=MPI_Allreduce((int*)&cnt,(int*)&res,1,MPI_INT,MPI_SUM,PETSC_COMM_WORLD);CHKERRQ(ierr);
00179 res*=5;
00180
00181 PetscPrintf(PETSC_COMM_WORLD,"sampling points: %i\n",res);
00182
00183
00184 for (i=0;i<cnt;i++) {
00185 for (j=0;j<NV;j++) {
00186 x2min[2]=PetscMin(x2min[2],vertxyz2[ND*gdata->elevert[NV*sliceele[i]+j]+2]);
00187 x2max[2]=PetscMax(x2max[2],vertxyz2[ND*gdata->elevert[NV*sliceele[i]+j]+2]);
00188 }
00189 }
00190
00191 ierr = PetscGlobalMin(&x2min[2],&xmin[2],PETSC_COMM_WORLD);CHKERRQ(ierr);
00192 x2min[2]=xmin[2];
00193 ierr = PetscGlobalMax(&x2max[2],&xmax[2],PETSC_COMM_WORLD);CHKERRQ(ierr);
00194 x2max[2]=xmax[2];
00195
00196 PetscPrintf(PETSC_COMM_WORLD,
00197 "sampling line ends in local coordinates (axes2):\n"
00198 "(%g,%g,%g) (%g,%g,%g)\n",
00199 x2min[0],x2min[1],x2min[2],
00200 x2max[0],x2max[1],x2max[2]
00201 );
00202
00203
00204 x2min[2] -= 1.0*(x2max[2]-x2min[2])/(res);
00205 x2max[2] += 2.0*(x2max[2]-x2min[2])/(res);
00206 PetscPrintf(PETSC_COMM_WORLD,
00207 "enlarged slice corners in local slice plane coordinates (axes2):\n"
00208 "(%g,%g,%g) (%g,%g,%g)\n",
00209 x2min[0],x2min[1],x2min[2],
00210 x2max[0],x2max[1],x2max[2]
00211 );
00212
00213 line_length=x2max[2]-x2min[2];
00214
00215
00216 AxesRotation(line_p,PETSC_NULL,axes2,1,1,x2min,xmin);
00217 AxesRotation(line_p,PETSC_NULL,axes2,1,1,x2max,xmax);
00218
00219 PetscPrintf(PETSC_COMM_WORLD,
00220 "sampling line ends in world coordinates (axes):\n"
00221 "(%g,%g,%g) (%g,%g,%g)\n",
00222 xmin[0],xmin[1],xmin[2],
00223 xmax[0],xmax[1],xmax[2]
00224 );
00225
00226 if (!rank) {
00227 PetscPrintf(PETSC_COMM_WORLD,"length of sampling line: %g\n",line_length);
00228
00229 char fmesh[256];
00230 FILE *fd;
00231 ierr = PetscSNPrintf(fmesh,255,"%s.%04i.datmsh",gdata->simname,gdata->inp);CHKERRQ(ierr);
00232 ierr = PetscFOpen(PETSC_COMM_WORLD,fmesh,"w",&fd);CHKERRQ(ierr);
00233 if (!fd) {
00234 SETERRQ1(PETSC_ERR_FILE_OPEN,"Cannot open %s\n",fmesh);
00235 }
00236 PetscPrintf(PETSC_COMM_WORLD,"Opening file %s\n",fmesh);
00237
00238 ierr = PetscFPrintf(PETSC_COMM_WORLD,fd,
00239 "#.%4s %14s %14s %14s %14s\n"
00240 "#.%4s %14s %14s %14s %14s\n",
00241 "..id","dist","x","y","z",
00242 "...-","-","-","-","-"
00243 );CHKERRQ(ierr);
00244
00245 for (i=0;i<res;i++) {
00246
00247 my_dcopy(ND,xmin,1,p,1);
00248 my_daxpy(ND,i*line_length/res,axes2+ND*2,1,p,1);
00249
00250 ierr = PetscFPrintf(PETSC_COMM_WORLD,fd,
00251 " %4d %14e %14e %14e %14e\n",
00252 i,i*line_length/res,
00253 p[0],p[1],p[2]
00254 );
00255 }
00256
00257 ierr = PetscFClose(PETSC_COMM_WORLD,fd);CHKERRQ(ierr);
00258 }
00259
00260
00261
00262 if (!rank) {
00263 ierr = MatCreateMPIAIJ(
00264 PETSC_COMM_WORLD,
00265 ND*res,ND*gdata->ln_vert,
00266 ND*res,ND*gdata->n_vert,
00267 NV+1,PETSC_NULL,NV+1,PETSC_NULL,
00268 &AIpolline
00269 );CHKERRQ(ierr);
00270 }
00271 else {
00272 ierr = MatCreateMPIAIJ(
00273 PETSC_COMM_WORLD,
00274 0,ND*gdata->ln_vert,
00275 ND*res,ND*gdata->n_vert,
00276 NV,PETSC_NULL,NV,PETSC_NULL,
00277 &AIpolline
00278 );CHKERRQ(ierr);
00279 }
00280 ierr = MatSetFromOptions(AIpolline);CHKERRQ(ierr);
00281
00282 if (cnt>0) {
00283 t_ele=sliceele[0];
00284 d=gdata->vertxyz+ND*gdata->elevert[NV*t_ele+3];
00285
00286 for (i=0;i<res;i++) {
00287
00288
00289 my_dcopy(ND,xmin,1,p,1);
00290 my_daxpy(ND,i*line_length/res,axes2+ND*2,1,p,1);
00291
00292
00293 for (k=-1;k<cnt;k++) {
00294
00295
00296
00297 if (k>=0) {
00298 t_ele=sliceele[k];
00299 d=gdata->vertxyz+ND*gdata->elevert[NV*t_ele+3];
00300 }
00301
00302
00303
00304
00305
00306
00307
00308
00309 my_dcopy(ND,p,1,rhs,1);
00310 my_daxpy(ND,-1.0,d,1,rhs,1);
00311
00312
00313 if (my_dnrm2(ND,rhs,1)>gdata->elenmax) continue;
00314
00315 if (
00316 barycent(
00317 p,
00318 gdata->vertxyz+ND*gdata->elevert[NV*t_ele+0],
00319 gdata->vertxyz+ND*gdata->elevert[NV*t_ele+1],
00320 gdata->vertxyz+ND*gdata->elevert[NV*t_ele+2],
00321 gdata->vertxyz+ND*gdata->elevert[NV*t_ele+3],
00322 rhs
00323 )
00324 ) {
00325 for (l=0;l<ND;l++) {
00326 matele=(1-rhs[0]-rhs[1]-rhs[2]);
00327 ierr = MatSetValue(AIpolline,
00328 ND*i+l,
00329 ND*gdata->elevert[NV*t_ele+3]+l,
00330 matele,ADD_VALUES);CHKERRQ(ierr);
00331
00332 matele=rhs[0];
00333 ierr = MatSetValue(AIpolline,
00334 ND*i+l,
00335 ND*gdata->elevert[NV*t_ele+0]+l,
00336 matele,ADD_VALUES);CHKERRQ(ierr);
00337
00338 matele=rhs[1];
00339 ierr = MatSetValue(AIpolline,
00340 ND*i+l,
00341 ND*gdata->elevert[NV*t_ele+1]+l,
00342 matele,ADD_VALUES);CHKERRQ(ierr);
00343
00344 matele=rhs[2];
00345 ierr = MatSetValue(AIpolline,
00346 ND*i+l,
00347 ND*gdata->elevert[NV*t_ele+2]+l,
00348 matele,ADD_VALUES);CHKERRQ(ierr);
00349 }
00350 break;
00351 }
00352 }
00353 }
00354 }
00355
00356
00357 PetscInfo(0,"AIpolline matrix assembly complete\n");
00358 ierr = MatAssemblyBegin(AIpolline,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
00359 ierr = MatAssemblyEnd(AIpolline,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
00360 ierr = PrintMatInfoAll(AIpolline);CHKERRQ(ierr);
00361
00362 ierr = VecCreate(PETSC_COMM_WORLD,&VIpolline);CHKERRQ(ierr);
00363 if (!rank) {
00364 ierr = VecSetSizes(VIpolline,ND*res,ND*res);CHKERRQ(ierr);
00365 }
00366 else {
00367 ierr = VecSetSizes(VIpolline,0,ND*res);CHKERRQ(ierr);
00368 }
00369 ierr = VecSetFromOptions(VIpolline);CHKERRQ(ierr);
00370
00371 ierr = PetscFree(sliceele);CHKERRQ(ierr);
00372 ierr = PetscFree(vertxyz2);CHKERRQ(ierr);
00373
00374
00375 MagparFunctionLogReturn(0);
00376 }
00377
00378
00379
00380 int WriteDat(GridData *gdata)
00381 {
00382 MagparFunctionInfoBegin;
00383
00384 int rank,size;
00385 ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRQ(ierr);
00386 ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRQ(ierr);
00387
00388 if (doinit) {
00389 ierr = WriteDatInit(gdata);CHKERRQ(ierr);
00390 doinit=0;
00391 }
00392 if (AIpolline==PETSC_NULL) {
00393 MagparFunctionInfoReturn(0);
00394 }
00395
00396 char fmesh[256];
00397 FILE *fd;
00398 ierr = PetscSNPrintf(fmesh,255,"%s.%04i.d",gdata->simname,gdata->inp);CHKERRQ(ierr);
00399 ierr = PetscFOpen(PETSC_COMM_WORLD,fmesh,"w",&fd);CHKERRQ(ierr);
00400 if (!rank && !fd) {
00401 SETERRQ1(PETSC_ERR_FILE_OPEN,"Cannot open %s\n",fmesh);
00402 }
00403 PetscPrintf(PETSC_COMM_WORLD,"Opening file %s\n",fmesh);
00404
00405 #ifdef ADDONS
00406
00407
00408 ierr = PetscFPrintf(PETSC_COMM_WORLD,fd,
00409 "#.%4s %14s %14s %14s\n"
00410 "#.%4s %14s %14s %14s\n",
00411 "..id","Htot_x","Htot_y","Htot_z",
00412 "....","(T)","(T)","(T)"
00413 );CHKERRQ(ierr);
00414
00415 ierr = MatMult(AIpolline,gdata->VHtot,VIpolline);CHKERRQ(ierr);
00416 ierr = VecScale(VIpolline,gdata->hscale);CHKERRQ(ierr);
00417
00418 #else
00419
00420 ierr = PetscFPrintf(PETSC_COMM_WORLD,fd,
00421 "#.%4s %14s %14s %14s\n"
00422 "#.%4s %14s %14s %14s %s\n",
00423 "..id","Mx","My","Mz",
00424 "...-","(1)","(1)","(1)",
00425 "(normalized magnetization, not proportional to Js)"
00426 );CHKERRQ(ierr);
00427
00428 ierr = MatMult(AIpolline,gdata->M,VIpolline);CHKERRQ(ierr);
00429 #endif
00430
00431 if (!rank) {
00432 int i,j;
00433 VecGetSize(VIpolline,(PetscInt*)&i);
00434 VecGetLocalSize(VIpolline,(PetscInt*)&j);
00435
00436 assert(i==j);
00437 PetscReal *ta_dat;
00438 VecGetArray(VIpolline,&ta_dat);
00439
00440 for (i=0; i<j/ND; i++) {
00441 ierr = PetscFPrintf(PETSC_COMM_WORLD,fd,
00442 " %4d %14e %14e %14e\n",
00443 i,
00444 ta_dat[ND*i+0],
00445 ta_dat[ND*i+1],
00446 ta_dat[ND*i+2]
00447 );
00448 }
00449 VecRestoreArray(VIpolline,&ta_dat);
00450 }
00451
00452 ierr = PetscFClose(PETSC_COMM_WORLD,fd);CHKERRQ(ierr);
00453
00454 MagparFunctionInfoReturn(0);
00455 }
00456