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