00001
00002
00003
00004
00005
00006 #ifndef AMROC_WENOSTATISTICS_H
00007 #define AMROC_WENOSTATISTICS_H
00008
00016 #include <iostream>
00017 #include <fstream>
00018 #include <string>
00019 #include <stdlib.h>
00020 #include "DAGH.h"
00021 #include "AMRBase.h"
00022 #include "Interpolation.h"
00023
00024 #include "StatParser/StatParser.h"
00025
00026
00027
00028 template <class VectorType, class InterpolationType, class OutputType, int dim>
00029 class WENOStatistics : public AMRBase<VectorType,dim>,
00030 private StatParser<typename InterpolationType::point_type,
00031 typename VectorType::InternalDataType,dim> {
00032 typedef AMRBase<VectorType,dim> base;
00033 typedef typename VectorType::InternalDataType DataType;
00034 typedef typename InterpolationType::point_type point_type;
00035 typedef StatParser<point_type,DataType,dim> stparser;
00036 typedef typename stparser::symrec symrec;
00037
00038
00039 using StatParser<point_type,DataType,dim>::group_table;
00040 using StatParser<point_type,DataType,dim>::putsym;
00041 using StatParser<point_type,DataType,dim>::unlinksym;
00042 using StatParser<point_type,DataType,dim>::vectorize;
00043 using StatParser<point_type,DataType,dim>::buildcoords;
00044 using StatParser<point_type,DataType,dim>::clearcoords;
00045
00046 public:
00047 typedef typename base::vec_grid_fct_type vec_grid_fct_type;
00048 typedef GridFunction<DataType,dim> grid_fct_type;
00049
00050 public:
00051 WENOStatistics(InterpolationType *int_ref, OutputType* out_ref) : base() {
00052 _Step = 1;
00053 _FileOutput = out_ref;
00054 _Interpolation = int_ref;
00055 _OutputFile = "stats";
00056 _DumpMeshes = 0;
00057 }
00058
00059 void Setup(GridHierarchy *gh, const int& ghosts, int* shape, double* geom) {
00060 base::SetupData(gh, ghosts);
00061
00062 int MaxLev;
00063 int nx, ny, nz;
00064 char s[36];
00065
00066 MaxLev = MaxLevel(base::GH());
00067
00068 nx = shape[0];
00069 putsym(TOK_INTEGER, STAT_NX)->vdata.ivalue = nx;
00070
00071 sprintf(s, "%s(%d)", STAT_NX, 1);
00072 putsym(TOK_INTEGER, s)->vdata.ivalue = nx;
00073 for ( int i = 0; i < MaxLev ; i++ ) {
00074 sprintf(s, "%s(%d)", STAT_NX, i+2);
00075 nx *= RefineFactor(base::GH(),i);
00076 putsym(TOK_INTEGER, s)->vdata.ivalue = nx;
00077 }
00078 putsym(TOK_INTEGER, STAT_NXF)->vdata.ivalue = nx;
00079 putsym(TOK_DOUBLE, STAT_XMIN)->vdata.dvalue = geom[0];
00080 putsym(TOK_DOUBLE, STAT_XMAX)->vdata.dvalue = geom[1];
00081
00082 if ( dim > 1 ) {
00083 ny = shape[1];
00084 putsym(TOK_INTEGER, STAT_NY)->vdata.ivalue = ny;
00085
00086 sprintf(s, "%s(%d)", STAT_NY, 1);
00087 putsym(TOK_INTEGER, s)->vdata.ivalue = ny;
00088 for ( int i = 0; i < MaxLev ; i++ ) {
00089 sprintf(s, "%s(%d)", STAT_NY, i+2);
00090 ny *= RefineFactor(base::GH(),i);;
00091 putsym(TOK_INTEGER, s)->vdata.ivalue = ny;
00092 }
00093 putsym(TOK_INTEGER, STAT_NYF)->vdata.ivalue = ny;
00094 putsym(TOK_DOUBLE, STAT_YMIN)->vdata.dvalue = geom[2];
00095 putsym(TOK_DOUBLE, STAT_YMAX)->vdata.dvalue = geom[3];
00096 }
00097
00098 if ( dim > 2 ) {
00099 nz = shape[2];
00100 putsym(TOK_INTEGER, STAT_NZ)->vdata.ivalue = nz;
00101
00102 sprintf(s, "%s(%d)", STAT_NZ, 1);
00103 putsym(TOK_INTEGER, s)->vdata.ivalue = nz;
00104 for ( int i = 0; i < MaxLev ; i++ ) {
00105 sprintf(s, "%s(%d)", STAT_NZ, i+2);
00106 nz *= RefineFactor(base::GH(),i);;
00107 putsym(TOK_INTEGER, s)->vdata.ivalue = nz;
00108 }
00109 putsym(TOK_INTEGER, STAT_NZF)->vdata.ivalue = nz;
00110 putsym(TOK_DOUBLE, STAT_ZMIN)->vdata.dvalue = geom[4];
00111 putsym(TOK_DOUBLE, STAT_ZMAX)->vdata.dvalue = geom[5];
00112 }
00113
00114
00115
00116
00117 symrec *ptr;
00118
00119 putsym(TOK_VARIABLE, "rho")->vdata.iptr = 1;
00120 putsym(TOK_VARIABLE, "r")->vdata.iptr = 1;
00121 (ptr = putsym(TOK_VARIABLE, "u"))->vdata.iptr = 2; ptr->src = STAT_DERIVED;
00122 if ( dim > 1 ) {
00123 (ptr = putsym(TOK_VARIABLE, "v"))->vdata.iptr = 3; ptr->src = STAT_DERIVED;
00124 }
00125 if ( dim > 2 ) {
00126 (ptr = putsym(TOK_VARIABLE, "w"))->vdata.iptr = 4; ptr->src = STAT_DERIVED;
00127 }
00128 putsym(TOK_VARIABLE, "E")->vdata.iptr = 5;
00129 for ( int i = 0; i < NVARS - 5 ; i++ ) {
00130 char name[12];
00131 sprintf(name, "Y%d", i+1);
00132 (ptr = putsym(TOK_VARIABLE, name))->vdata.iptr =
00133 ((dim == 1 ? 7 : (dim == 2 ? 8 : 9))+i); ptr->src = STAT_DERIVED;
00134 }
00135 if ( dim == 3 ) {
00136 (ptr = putsym(TOK_VARIABLE, "sgske"))->vdata.iptr = (NVARS+5); ptr->src = STAT_DERIVED;
00137 }
00138 putsym(TOK_VARIABLE, "T")->vdata.iptr = NVARS+1;
00139 (ptr = putsym(TOK_VARIABLE, "p"))->vdata.iptr =
00140 (dim == 1 ? 5 : (dim == 2 ? 6 : 7)); ptr->src = STAT_DERIVED;
00141 (ptr = putsym(TOK_VARIABLE, "gamma"))->vdata.iptr =
00142 (dim == 1 ? 6 : (dim == 2 ? 7 : 8)); ptr->src = STAT_DERIVED;
00143 putsym(TOK_VARINARRAY, "q")->vdata.iptr = 0;
00144 (ptr = putsym(TOK_VARINARRAY, "qout"))->vdata.iptr = 0; ptr->src = STAT_DERIVED;
00145
00146
00147 if ( _DefFile.length() ) {
00148 stparser::parse(_DefFile.c_str());
00149
00150 for ( typename stparser::grpvec::const_iterator g = group_table.begin() ;
00151 g != group_table.end(); g++)
00152 if ( (*g)->step == 0 ) (*g)->step = _Step;
00153
00154
00155 int me = MY_PROC;
00156 if (me == VizServer) {
00157 if ( _DumpMeshes ) {
00158 int ig = 1;
00159 for ( typename stparser::grpvec::const_iterator g = group_table.begin() ;
00160 g != group_table.end(); g++, ig++) {
00161
00162 int ip = 1;
00163 for ( typename stparser::provec::const_iterator p = (*g)->probe_table.begin() ;
00164 p != (*g)->probe_table.end() ; p ++, ip++ )
00165 {
00166 if ( buildcoords(*p) != RETURN_OK ) return;
00167
00168 int i, j;
00169 FILE * fp;
00170 char name[128];
00171
00172 sprintf(name, "group%d-probe%d.dat", ig, ip);
00173 fp = fopen(name, "w");
00174 if ( (*p)->type == TOK_THREAD ) {
00175 for ( i = 0 ; i < (*p)->options.thread.na ; i++ ) {
00176 fprintf(fp, "%12.8lf %12.8lf %12.8lf\n", (*p)->coords[i](0),
00177 (*p)->coords[i](1),
00178 (*p)->coords[i](2));
00179 }
00180 } else {
00181 int na, nb;
00182 if ( (*p)->type == TOK_SURFACE ) {
00183 na = (*p)->options.surface.na;
00184 nb = (*p)->options.surface.nb;
00185 } else if ((*p)->type == TOK_PLANES ) {
00186 na = (*p)->options.plane.na;
00187 nb = (*p)->options.plane.nb;
00188 }
00189 fprintf(fp, "VARIABLES = \"X\", \"Y\", \"Z\"\n");
00190 fprintf(fp, "ZONE I=%d, J=%d, F=POINT\n", na, nb);
00191 for ( j = 0 ; j < nb ; j++ )
00192 for ( i = 0 ; i < na ; i++ ) {
00193 int l = i+j*na;
00194 fprintf(fp, "%12.8lf %12.8lf %12.8lf\n", (*p)->coords[l](0),
00195 (*p)->coords[l](1), (*p)->coords[l](2));
00196 }
00197 }
00198 fclose(fp);
00199
00200 clearcoords(*p);
00201 }
00202 }
00203 }
00204 }
00205 }
00206 }
00207
00208 virtual void register_at(ControlDevice& Ctrl) { register_at(Ctrl, ""); }
00209 virtual void register_at(ControlDevice& Ctrl, const std::string& prefix) {
00210 LocCtrl = Ctrl.getSubDevice(prefix+"Statistics");
00211
00212 RegisterAt(LocCtrl,"DefinitionFile",_DefFile);
00213 RegisterAt(LocCtrl,"OutputFile",_OutputFile);
00214 RegisterAt(LocCtrl,"Step",_Step);
00215 RegisterAt(LocCtrl,"WriteMeshes", _DumpMeshes);
00216 }
00217
00218 ~WENOStatistics() {
00219 _FileOutput = NULL;
00220 _Interpolation = NULL;
00221 }
00222
00223 public:
00224
00225 void Evaluate(double * t, vec_grid_fct_type& U, grid_fct_type& Work)
00226 {
00227 int me ;
00228 int Npoints;
00229
00230 char fname[256];
00231 std::ofstream outfile;
00232 std::ostream* out;
00233
00234 if ( group_table.size() == 0 ) return;
00235
00236
00237 int TimeC = StepsTaken(base::GH(),0);
00238
00239
00240 int gcheck = 0;
00241 for ( typename stparser::grpvec::const_iterator g = group_table.begin() ;
00242 g != group_table.end(); g++) {
00243 if ( TimeC%((*g)->step) == 0 ) {
00244 gcheck = 1;
00245 break;
00246 }
00247 }
00248
00249
00250 if ( gcheck == 0 ) return;
00251
00252 me = MY_PROC;
00253
00254
00255 if (me == VizServer) {
00256 int time = CurrentTime(base::GH(), 0);
00257 std::sprintf(fname,"%s%d.dat",_OutputFile.c_str(),time);
00258 outfile.open(fname, std::ios::out );
00259 out = new std::ostream(outfile.rdbuf());
00260
00261 *out << "# time " << t[0] << std::endl;
00262 }
00263
00264
00265 int ig = 1, iptr;
00266 for ( typename stparser::grpvec::const_iterator g = group_table.begin() ;
00267 g != group_table.end(); g++, ig++ ) {
00268
00269 if ( TimeC%((*g)->step) == 0 ) {
00270
00271 if (me == VizServer)
00272 *out << "# group " << ig << std::endl;
00273
00274 Npoints = 0;
00275 for ( typename stparser::provec::const_iterator p = (*g)->probe_table.begin() ;
00276 p != (*g)->probe_table.end() ; p ++) {
00277 switch ( (*p)->type ) {
00278 case TOK_THREAD:
00279 Npoints += (*p)->options.thread.na;
00280 break;
00281 case TOK_SURFACE:
00282 Npoints += (*p)->options.surface.na*(*p)->options.surface.nb;
00283 break;
00284 case TOK_PLANES:
00285 Npoints += (*p)->options.plane.na*(*p)->options.plane.nb;
00286 break;
00287 case TOK_VOLUME:
00288
00289 break;
00290 }
00291 }
00292
00293 point_type* CList = new point_type[Npoints];
00294 int cli=0;
00295
00296 for ( typename stparser::provec::const_iterator p = (*g)->probe_table.begin() ;
00297 p != (*g)->probe_table.end() ; p ++) {
00298 if ( buildcoords(*p) != RETURN_OK ) return;
00299 int Lpoints;
00300 switch ( (*p)->type ) {
00301 case TOK_THREAD:
00302 Lpoints = (*p)->options.thread.na;
00303 break;
00304 case TOK_SURFACE:
00305 Lpoints = (*p)->options.surface.na*(*p)->options.surface.nb;
00306 break;
00307 case TOK_PLANES:
00308 Lpoints = (*p)->options.plane.na*(*p)->options.plane.nb;
00309 break;
00310 case TOK_VOLUME:
00311
00312 break;
00313 }
00314 for (register int li=0; li<Lpoints; li++, cli++)
00315 CList[cli] = (*p)->coords[li];
00316 if (me == VizServer)
00317 (*p)->out = new std::stringstream(std::stringstream::in | std::stringstream::out);
00318 }
00319
00320 bool *lused = new bool[Npoints];
00321 int Level = 0;
00322 int Time = CurrentTime(base::GH(),Level);
00323 int NpUsed =_Interpolation->LocalPoints(Work,Time,Level,Npoints,CList,lused);
00324 point_type* CUList = new point_type[NpUsed];
00325 cli = 0;
00326 for (register int i=0; i<Npoints; i++)
00327 if (lused[i]) CUList[cli++] = CList[i];
00328 delete [] CList;
00329
00330 symrec *ptr, *ptr_ref;
00331
00332
00333 ptr = stparser::clonestack((*g)->sym_table);
00334
00335
00336 while ( ptr ) {
00337 #ifdef DEBUG_PRINT_STATISTICS
00338 (comm_service::log() << stparser::token(ptr->type) << "\n").flush();
00339 #endif
00340
00341 vectorize(ptr, NpUsed);
00342
00343
00344
00345
00346 switch( ptr->type ) {
00347 case TOK_OP_PLUS:
00348 vectorize(ptr->prev, NpUsed);
00349 vectorize(ptr->prev->prev, NpUsed);
00350 for ( register int l = 0 ; l < NpUsed ; l++ )
00351 ptr->pdata[l] = ptr->prev->prev->pdata[l] + ptr->prev->pdata[l];
00352 unlinksym(ptr->prev->prev);
00353 unlinksym(ptr->prev);
00354 break;
00355 case TOK_OP_MINUS:
00356 vectorize(ptr->prev, NpUsed);
00357 vectorize(ptr->prev->prev, NpUsed);
00358 for ( register int l = 0 ; l < NpUsed ; l++ )
00359 ptr->pdata[l] = ptr->prev->prev->pdata[l] - ptr->prev->pdata[l];
00360 unlinksym(ptr->prev->prev);
00361 unlinksym(ptr->prev);
00362 break;
00363 case TOK_OP_PROD:
00364 vectorize(ptr->prev, NpUsed);
00365 vectorize(ptr->prev->prev, NpUsed);
00366 for ( register int l = 0 ; l < NpUsed ; l++ )
00367 ptr->pdata[l] = ptr->prev->prev->pdata[l] * ptr->prev->pdata[l];
00368 unlinksym(ptr->prev->prev);
00369 unlinksym(ptr->prev);
00370 break;
00371 case TOK_OP_DIV:
00372 vectorize(ptr->prev, NpUsed);
00373 vectorize(ptr->prev->prev, NpUsed);
00374 for ( register int l = 0 ; l < NpUsed ; l++ )
00375 ptr->pdata[l] = ptr->prev->prev->pdata[l] / ptr->prev->pdata[l];
00376 unlinksym(ptr->prev->prev);
00377 unlinksym(ptr->prev);
00378 break;
00379 case TOK_OP_NEG:
00380 vectorize(ptr->prev, NpUsed);
00381 for ( register int l = 0 ; l < NpUsed ; l++ )
00382 ptr->pdata[l] = - ptr->prev->pdata[l];
00383 unlinksym(ptr->prev);
00384 break;
00385 case TOK_OP_POWER:
00386 vectorize(ptr->prev, NpUsed);
00387 vectorize(ptr->prev->prev, NpUsed);
00388 for ( register int l = 0 ; l < NpUsed ; l++ )
00389 ptr->pdata[l] = pow(ptr->prev->prev->pdata[l], ptr->prev->pdata[l]);
00390 unlinksym(ptr->prev->prev);
00391 unlinksym(ptr->prev);
00392 break;
00393 case TOK_POINTER:
00394 ptr_ref = ptr->ref;
00395 switch(ptr_ref->type) {
00396 case TOK_FUNCTION1:
00397 vectorize(ptr->prev, NpUsed);
00398 for ( register int l = 0 ; l < NpUsed ; l++ )
00399 ptr->pdata[l] = ((func_t1)ptr_ref->vdata.ptr)(ptr->prev->pdata[l]);
00400 unlinksym(ptr->prev);
00401 break;
00402 case TOK_FUNCTION2:
00403 vectorize(ptr->prev, NpUsed);
00404 vectorize(ptr->prev->prev, NpUsed);
00405 for ( register int l = 0 ; l < NpUsed ; l++ )
00406 ptr->pdata[l] = ((func_t2)ptr_ref->vdata.ptr)(ptr->prev->prev->pdata[l],
00407 ptr->prev->pdata[l]);
00408 unlinksym(ptr->prev->prev);
00409 unlinksym(ptr->prev);
00410 break;
00411 case TOK_VARINARRAY:
00412 case TOK_VARIABLE:
00413
00414
00415
00416
00417 if ( ptr_ref->type == TOK_VARINARRAY )
00418 iptr = ptr->vdata.iptr;
00419 else
00420 iptr = ptr_ref->vdata.iptr;
00421
00422
00423 if ( ptr_ref->src == STAT_DERIVED ) {
00424
00425 for (register int l=0; l<=FineLevel(base::GH()); l++) {
00426 int Time = CurrentTime(base::GH(),l);
00427
00428 _FileOutput->Transform(U, Work, Time, l, iptr, t[l]);
00429 }
00430 #ifdef DEBUG_PRINT_STATISTICS
00431 (comm_service::log() << " convert output\n" ).flush();
00432 #endif
00433 } else {
00434 int _dim = dim;
00435
00436 if ( _dim == 2 ) {
00437 for ( int Level = FineLevel(base::GH()) ; Level >= 0 ; Level --) {
00438 int Time = CurrentTime(base::GH(),Level);
00439 forall (U,Time,Level,c)
00440 Coords ss(U(Time,Level,c).bbox().stepsize());
00441 BBox bbi(U.boundingbbox(Time,Level,c));
00442 BeginFastIndex2(U, U(Time,Level,c).bbox(),
00443 U(Time,Level,c).data(), VectorType);
00444 BeginFastIndex2(Work, Work(Time,Level,c).bbox(),
00445 Work(Time,Level,c).data(), DataType);
00446 for_2 (i, j, bbi, ss)
00447 FastIndex2(Work,i,j) = FastIndex2(U,i,j)(iptr-1);
00448 end_for
00449 EndFastIndex2(U);
00450 EndFastIndex2(Work);
00451 end_forall
00452 }
00453 } else if ( _dim == 3 ) {
00454 for ( int Level = FineLevel(base::GH()) ; Level >= 0 ; Level --) {
00455 int Time = CurrentTime(base::GH(),Level);
00456 forall (U,Time,Level,c)
00457 Coords ss(U(Time,Level,c).bbox().stepsize());
00458 BBox bbi(U.boundingbbox(Time,Level,c));
00459 BeginFastIndex3(U, U(Time,Level,c).bbox(),
00460 U(Time,Level,c).data(), VectorType);
00461 BeginFastIndex3(Work, Work(Time,Level,c).bbox(),
00462 Work(Time,Level,c).data(), DataType);
00463 for_3 (i, j, k, bbi, ss)
00464 FastIndex3(Work,i,j,k) = FastIndex3(U,i,j,k)(iptr-1);
00465 end_for
00466 EndFastIndex3(U);
00467 EndFastIndex3(Work);
00468 end_forall
00469 }
00470 }
00471 #ifdef DEBUG_PRINT_STATISTICS
00472 (comm_service::log() << " convert array\n" ).flush();
00473 #endif
00474 }
00475
00476 _Interpolation->PointsValues(Work,t[0],NpUsed,CUList,ptr->pdata);
00477 #ifdef DEBUG_PRINT_STATISTICS
00478 (comm_service::log() << "after\n" ).flush();
00479 #endif
00480 break;
00481 case TOK_COORDINATE:
00482 {
00483 int iptr = (ptr_ref->vdata.iptr)-1;
00484 if ( iptr >= dim ) {
00485 staterror("invalid coordinate in this stack");
00486 return;
00487 }
00488 for ( register int l = 0 ; l < NpUsed ; l++ )
00489 ptr->pdata[l] = CUList[l](iptr);
00490 }
00491 break;
00492 default:
00493 staterror("there should not be any other pointer type left in this stack");
00494 return;
00495 }
00496
00497 break;
00498
00499 case TOK_SEPARATOR:
00500
00501 {
00502 int na, nb, avea, aveb;
00503 vectorize(ptr->prev, NpUsed);
00504
00505 cli=0;
00506 int ip = 1;
00507 int ccli=0;
00508 for ( typename stparser::provec::const_iterator p = (*g)->probe_table.begin() ;
00509 p != (*g)->probe_table.end() ; p ++, ip++ ) {
00510
00511 int Lpoints;
00512 switch((*p)->type) {
00513 case TOK_THREAD:
00514 Lpoints = (*p)->options.thread.na;
00515 if ( (*p)->options.thread.ave == STAT_AVERAGED ) {
00516 double sum = 0.0, rsum;
00517 for ( register int ip = 0 ; ip < Lpoints ; ip ++ )
00518 if (lused[ip+cli])
00519 sum += ptr->prev->pdata[ccli++];
00520 MPI_Reduce(&sum, &rsum, 1, MPI_DOUBLE, MPI_SUM, VizServer, comm_service::comm());
00521 if (me == VizServer)
00522 *(*p)->out << (rsum/Lpoints) << std::endl;
00523 } else {
00524 double *adata = new double[Lpoints];
00525 for ( register int ip = 0 ; ip < Lpoints ; ip ++ )
00526 if (lused[ip+cli])
00527 adata[ip] = ptr->prev->pdata[ccli++];
00528 else
00529 adata[ip] = _Interpolation->ErrorValue();
00530 _Interpolation->ArrayCombine(VizServer, Lpoints, adata);
00531 if (me == VizServer) {
00532 for ( register int ip = 0 ; ip < Lpoints ; ip ++ )
00533 *(*p)->out << adata[ip] << " ";
00534 *(*p)->out << std::endl;
00535 }
00536 delete [] adata;
00537 }
00538 break;
00539 case TOK_SURFACE:
00540 Lpoints = (*p)->options.surface.na*(*p)->options.surface.nb;
00541 avea = (*p)->options.surface.avea;
00542 aveb = (*p)->options.surface.aveb;
00543 na = (*p)->options.surface.na;
00544 nb = (*p)->options.surface.nb;
00545 case TOK_PLANES:
00546 if ( (*p)->type == TOK_PLANES ) {
00547 Lpoints = (*p)->options.plane.na*(*p)->options.plane.nb;
00548 avea = STAT_AVERAGED;
00549 aveb = STAT_AVERAGED;
00550 na = (*p)->options.plane.na;
00551 nb = (*p)->options.plane.nb;
00552 }
00553
00554 if ( avea == STAT_AVERAGED && aveb == STAT_AVERAGED ) {
00555 double sum = 0.0, rsum;
00556 for ( register int ip = 0 ; ip < Lpoints ; ip ++ )
00557 if (lused[ip+cli])
00558 sum += ptr->prev->pdata[ccli++];
00559 MPI_Reduce(&sum, &rsum, 1, MPI_DOUBLE, MPI_SUM, VizServer, comm_service::comm());
00560 if (me == VizServer)
00561 *(*p)->out << (rsum/Lpoints) << std::endl;
00562 } else {
00563 double *adata = new double[Lpoints];
00564 for ( register int j = 0 ; j < nb ; j++ )
00565 for ( register int i = 0 ; i < na ; i++ )
00566 if (lused[i+j*na+cli])
00567 adata[i+j*na] = ptr->prev->pdata[ccli++];
00568 else
00569 adata[i+j*na] = _Interpolation->ErrorValue();
00570
00571 if ( avea == STAT_RAW && aveb == STAT_AVERAGED ) {
00572 for ( register int i = 0 ; i < na ; i++ ) {
00573 double sum=0.0, rsum;
00574 for ( register int j = 0 ; j < nb ; j++ )
00575 if (adata[i+j*na]!=_Interpolation->ErrorValue())
00576 sum += adata[i+j*na];
00577 MPI_Reduce(&sum, &rsum, 1, MPI_DOUBLE, MPI_SUM, VizServer, comm_service::comm());
00578 if (me == VizServer)
00579 *(*p)->out << (rsum/nb) << " ";
00580 }
00581 if (me == VizServer)
00582 *(*p)->out << std::endl;
00583 } else if ( avea == STAT_AVERAGED && aveb == STAT_RAW ) {
00584 for ( register int j = 0 ; j < nb ; j++ ) {
00585 double sum=0.0, rsum;
00586 for ( register int i = 0 ; i < na ; i++ )
00587 if (adata[i+j*na]!=_Interpolation->ErrorValue())
00588 sum += adata[i+j*na];
00589 MPI_Reduce(&sum, &rsum, 1, MPI_DOUBLE, MPI_SUM, VizServer, comm_service::comm());
00590 if (me == VizServer)
00591 *(*p)->out << (rsum/na) << " ";
00592 }
00593 if (me == VizServer)
00594 *(*p)->out << std::endl;
00595 } else {
00596 _Interpolation->ArrayCombine(VizServer, Lpoints, adata);
00597 if (me == VizServer) {
00598 for ( register int ip = 0 ; ip < Lpoints ; ip ++ )
00599 *(*p)->out << adata[ip] << " ";
00600 *(*p)->out << std::endl;
00601 }
00602 }
00603 delete [] adata;
00604 }
00605 break;
00606 case TOK_VOLUME:
00607 break;
00608 }
00609 cli += Lpoints;
00610 clearcoords(*p);
00611 }
00612 unlinksym(ptr->prev);
00613 }
00614 #ifndef DAGH_NO_MPI
00615 MPI_Barrier(comm_service::comm());
00616 #endif
00617 break;
00618 case TOK_FUNCTION1:
00619 case TOK_FUNCTION2:
00620 case TOK_VARINARRAY:
00621 case TOK_VARIABLE:
00622 case TOK_COORDINATE:
00623 staterror("there should not be any pointer type tokens"
00624 " left in the stack");
00625 return;
00626 default:
00627 if ( ptr->type != TOK_DOUBLE ) {
00628 staterror("unhandled stack token");
00629 return;
00630 }
00631 break;
00632 }
00633
00634 if ( ptr->type == TOK_SEPARATOR ) {
00635 ptr_ref = ptr->next;
00636 unlinksym(ptr);
00637 ptr = ptr_ref;
00638 } else {
00639 ptr->type = TOK_DOUBLE;
00640 ptr = ptr->next;
00641 }
00642
00643 #ifdef DEBUG_PRINT_STATISTICS
00644 (comm_service::log() << " end of token switch\n" ).flush();
00645 #endif
00646 }
00647 delete [] lused;
00648 delete [] CUList;
00649
00650 if (me == VizServer) {
00651 int ip = 1;
00652 for ( typename stparser::provec::const_iterator p = (*g)->probe_table.begin() ;
00653 p != (*g)->probe_table.end() ; p ++, ip++) {
00654 *out << "# probe " << stparser::token((*p)->type) << " " << ip << std::endl;
00655 ( *out << (*p)->out->str() ).flush();;
00656 delete (*p)->out;
00657 }
00658 }
00659 }
00660 }
00661
00662
00663 if (me == VizServer) {
00664 outfile.close();
00665 delete out;
00666 }
00667 }
00668
00669 private:
00670 int _DumpMeshes;
00671 int _Step;
00672 ControlDevice LocCtrl;
00673 std::string _OutputFile;
00674 std::string _DefFile;
00675
00676 OutputType* _FileOutput;
00677 InterpolationType* _Interpolation;
00678 };
00679
00680 #endif