amroc/weno/WENOStatistics.h

Go to the documentation of this file.
00001 // -*- C++ -*-
00002 
00003 // Copyright (C) California Institute of Technology
00004 // Carlos Pantano
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 //#define DEBUG_PRINT_STATISTICS
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   // bring members from templated base class into scope
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     // add level dependent point number
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         // add level dependent point number
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         // add level dependent point number
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     /* simulation variables (dimension dependent) */
00115     // by convention I take positive iptr's as conservative
00116     // variables and negative iptr's as primitive variables.
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     // Process the rack definition file
00147     if ( _DefFile.length() ) { 
00148       stparser::parse(_DefFile.c_str());
00149       // if any group did not specify step sampling rate, then set it to default common value
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       // write probe information: coordinates
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             // for each probe
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       //int TimeC = CurrentTime(base::GH(),0);
00237       int TimeC = StepsTaken(base::GH(),0);
00238 
00239       // check step sampling rate
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       //if ( TimeC%_Step != 0 ) return;
00250       if ( gcheck == 0 ) return;
00251 
00252       me = MY_PROC; 
00253 
00254       // open output file
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         // write time header
00261         *out << "# time " << t[0] << std::endl;
00262       }
00263       
00264       // for each group
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           // for each probe
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               /* not implemented yet */
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               /* not implemented yet */
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           /* first, make a copy of the stack */
00333           ptr = stparser::clonestack((*g)->sym_table);
00334   
00335           /* traverse the stack and evaluate */
00336           while ( ptr ) {
00337 #ifdef DEBUG_PRINT_STATISTICS
00338             (comm_service::log() << stparser::token(ptr->type) << "\n").flush();
00339 #endif
00340             // vectorize the data array that will be used for the result
00341             vectorize(ptr, NpUsed);
00342             // the arguments of the operator do not need to be vectorized
00343             // for all processors since they are not used. Only VizServer
00344             // performs operations with them.
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                 // for variables in one of the AMR arrays, the pointer is located in the
00415                 // current token, not the defered one (ref). For variables in regular
00416                 // variables, the pointer is stored in the pointed variable.
00417                 if ( ptr_ref->type == TOK_VARINARRAY ) 
00418                   iptr = ptr->vdata.iptr;
00419                 else 
00420                   iptr = ptr_ref->vdata.iptr;
00421                 
00422                 // process the variables here
00423                 if ( ptr_ref->src == STAT_DERIVED ) {
00424                   // output variable type case
00425                   for (register int l=0; l<=FineLevel(base::GH()); l++) {
00426                     int Time = CurrentTime(base::GH(),l);
00427                     // if derived variable (call Transform)
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                   // conservative vector of state variable case
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                 // note that we use data so the pointer token will be converted.
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               // output
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                   // write probe data
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             // need special code to handle the separator token, to flush the queue
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       // close output file
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