amroc/converters/VisualGridData.h

Go to the documentation of this file.

// -*- C++ -*-

// Copyright (C) 2002 Ralf Deiterding
// Brandenburgische Universitaet Cottbus
//
// Copyright (C) 2003-2007 California Institute of Technology
// Ralf Deiterding, ralf@cacr.caltech.edu

#ifndef AMROC_VISUALGRIDDATA_H
#define AMROC_VISUALGRIDDATA_H

#include "GridData3.h" 
#include "Evaluator.h"

#include <iostream>
#include <vector>
#include <cfloat>
#include <cmath>

template <class VectorType> class VisualGrid;

template <class VectorType>
class VisualGridData {
  typedef typename VectorType::InternalDataType DataType;
  typedef VisualGridData<VectorType> self;

public:
  typedef VisualGrid<VectorType> grid_type;
  typedef Evaluator<VectorType> evaluator_type;
  typedef GridData<DataType,3> data_block_type;
  typedef GridData<VectorType,3> vector_block_type;

  typedef std::vector<self *> gblock_list_type;
  typedef typename gblock_list_type::iterator gblock_list_iterator;
  typedef std::vector<float> point_list_type;
  typedef typename point_list_type::iterator point_list_iterator;

  enum direction {x_dir, y_dir, z_dir};

public:
  VisualGridData(grid_type& gr, data_block_type& db, BBox &bb, int comp, evaluator_type& ev,
                 int lev, int proc, float* pdx, float* pgeom) : 
    _Grid(gr), _Evaluator(ev), _Level(lev), _Processor(proc), dx(pdx), geom(pgeom) {

    _data_block = new vector_block_type(bb);    
    BeginFastIndex3(source, db.bbox(), db.data(), DataType);
    BeginFastIndex3(target, _data_block->bbox(), _data_block->data(),
                    VectorType);
    if (comp>=0) {
      for_3 (i, j, k, bb, bb.stepsize())
        FastIndex3(target,i,j,k)(comp) = FastIndex3(source,i,j,k); 
      end_for
    }
    else {
      for_3 (i, j, k, bb, bb.stepsize())
        FastIndex3(target,i,j,k) = FastIndex3(source,i,j,k); 
      end_for
    }
    EndFastIndex3(source);
    EndFastIndex3(target);
    DimUsed = gr.Dimensions();
  }

  VisualGridData(grid_type& gr, vector_block_type& db, BBox &bb, evaluator_type& ev,
                 int lev, int proc, float* pdx, float* pgeom) : 
    _Grid(gr), _Evaluator(ev), _Level(lev), _Processor(proc), dx(pdx), geom(pgeom) {

    _data_block = new vector_block_type(bb);
    BeginFastIndex3(source, db.bbox(), db.data(), VectorType);
    BeginFastIndex3(target, _data_block->bbox(), _data_block->data(),
                    VectorType);
    for_3 (i, j, k, bb, bb.stepsize())
      FastIndex3(target,i,j,k) = FastIndex3(source,i,j,k); 
    end_for
    EndFastIndex3(source);
    EndFastIndex3(target);
    DimUsed = gr.Dimensions();
  }

  virtual ~VisualGridData() 
  { delete _data_block; }

public:
  virtual void SetBlocks(int blocks[],int& offset,int& nOffset) {
    std::cerr << ".";
    _nodeOffset = nOffset;
    blocks[offset]   = DB().extents(0)+DimUsed(0);
    blocks[offset+1] = DB().extents(1)+DimUsed(1);
    blocks[offset+2] = DB().extents(2)+DimUsed(2);
    nOffset += blocks[offset]*blocks[offset+1]*blocks[offset+2];   
    offset += 3;
  }

  virtual void SetGridCells(float xyz[][3],int& offset){
    std::cerr << ".";
    int i, j, k,orii,orij,orik,ssi,ssj,ssk,ijk,idx;
    int exti,extj,extk;
    exti=DB().extents(0);
    extj=DB().extents(1);
    extk=DB().extents(2);
    orii=DB().lower(0);
    orij=DB().lower(1);
    orik=DB().lower(2);
    ssi=DB().stepsize(0);
    ssj=DB().stepsize(1);
    ssk=DB().stepsize(2);
    
    idx = 0;
    for (k = 0; k < extk; k++){
      for (j = 0; j < extj; j++){
        for (i = 0; i < exti; i++){
          
          ijk = k*extj*exti + j*exti + i;
          
          xyz[ijk+offset][0] = (orii+((float)i+0.5)*ssi)*dx[0]+geom[0];
          xyz[ijk+offset][1] = (orij+((float)j+0.5)*ssj)*dx[1]+geom[2];
          xyz[ijk+offset][2] = (orik+((float)k+0.5)*ssk)*dx[2]+geom[4];
          idx++;
        }
      }
    }
    offset+=idx;
  }
  virtual void SetGridNodes(float xyz[][3],int& offset) {
    std::cerr << ".";
    int i, j, k,orii,orij,orik,ssi,ssj,ssk,ijk,idx;
    int exti,extj,extk;
    exti=DB().extents(0)+DimUsed(0);
    extj=DB().extents(1)+DimUsed(1);
    extk=DB().extents(2)+DimUsed(2);
    orii=DB().lower(0);
    orij=DB().lower(1);
    orik=DB().lower(2);
    ssi=DB().stepsize(0);
    ssj=DB().stepsize(1);
    ssk=DB().stepsize(2);
    
    idx = 0;
    for (k = 0; k < extk; k++){
      for (j = 0; j < extj; j++){
        for (i = 0; i < exti; i++){
          
          ijk = k*extj*exti + j*exti + i;
          
          xyz[ijk+offset][0] = (float)(orii+i*ssi)*dx[0]+geom[0];
          xyz[ijk+offset][1] = (float)(orij+j*ssj)*dx[1]+geom[2];
          xyz[ijk+offset][2] = (float)(orik+k*ssk)*dx[2]+geom[4];
          idx++;
        }
      }
    }
    offset+=idx;
  }
  
  virtual void SetGridConns(int con[][8],int& offset) {
    std::cerr << ".";
    int i,j,k,ijk,idx;
    int exti,extj,extk;
    exti=DB().extents(0)+DimUsed(0);
    extj=DB().extents(1)+DimUsed(1);
    extk=DB().extents(2)+DimUsed(2);
    
    idx = 0;
    for (k = 0; k < extk-DimUsed(2); k++){
      for (j = 0; j < extj-DimUsed(1); j++){
        for (i = 0; i < exti-DimUsed(0); i++){
          
          ijk = k*(extj-DimUsed(1))*(exti-DimUsed(0)) +
                j*(exti-DimUsed(0)) + i;
          
          con[ijk+offset][0] =  k*extj*exti +
                                j*exti +
                                i + _nodeOffset;
          con[ijk+offset][1] =  k*extj*exti +
                                j*exti +
                                (i+DimUsed(0)) + _nodeOffset;
          con[ijk+offset][2] =  k*extj*exti +
                                (j+DimUsed(1))*exti +
                                i + _nodeOffset; 
          con[ijk+offset][3] =  k*extj*exti +
                                (j+DimUsed(1))*exti +
                                (i+DimUsed(0)) + _nodeOffset; 
          con[ijk+offset][4] = (k+DimUsed(2))*extj*exti +
                               j*exti +
                               i + _nodeOffset;
          con[ijk+offset][5] = (k+DimUsed(2))*extj*exti +
                               j*exti +
                               (i+DimUsed(0)) + _nodeOffset; 
          con[ijk+offset][6] = (k+DimUsed(2))*extj*exti +
                               (j+DimUsed(1))*exti +
                               i + _nodeOffset; 
          con[ijk+offset][7] = (k+DimUsed(2))*extj*exti +
                               (j+DimUsed(1))*exti +
                               (i+DimUsed(0)) + _nodeOffset; 
          
          idx++;
        }
      }
    }
    offset+=idx;
  }
  
  virtual void SetScalCells(int *key,float v[],int& offset,bool* CompRead) {
    std::cerr << ".";
    if (*key==1 || *key==2) {
      int idx=0;
      for_3 (i,j,k,DB().bbox(),DB().bbox().stepsize())
        if (*key==1) {
          v[idx+offset]= Processor() + 1.0;
          TheEvaluator().mvals[0] = Min(TheEvaluator().mvals[0],v[idx+offset]);
          TheEvaluator().mvals[1] = Max(TheEvaluator().mvals[1],v[idx+offset]);
        }
        else {
          v[idx+offset]= Level() + 1.0;   
          TheEvaluator().mvals[2] = Min(TheEvaluator().mvals[2],v[idx+offset]);
          TheEvaluator().mvals[3] = Max(TheEvaluator().mvals[3],v[idx+offset]);
        }
        idx++;        
      end_for
      offset+=idx;
    }
    else
      TheEvaluator().SetScalCells(*key,v,offset,DB(),dx,CompRead);
  }  

  virtual void SetScalNodes(int *key,float v[],int& offset,bool* CompRead) {
    std::cerr << ".";
    if (*key==1) {
      int idx=0;
      BBox lbbox(grow(DB().bbox(),DimUsed));
      data_block_type dbwork(lbbox);
      dbwork = -1.0;
      BeginFastIndex3(tgt, lbbox, dbwork.data(), DataType);
      
      for (gblock_list_iterator bit=_Grid.BlockList().begin();
           bit!=_Grid.BlockList().end(); bit++) 
        if ((*bit)->Processor() == Processor()) {
          BBox bbcur(growupper((*bit)->DB().bbox(),1));
          bbcur.setstepsize(lbbox.stepsize());
          bbcur.growupper(-1);
          BBox bbwork(lbbox * bbcur);
          if (bbwork.empty())
            continue;
          
          bbwork.lower() = (bbwork.lower()/bbwork.stepsize())*
                           bbwork.stepsize();
          bbwork.upper() = (bbwork.upper()/bbwork.stepsize())*
                           bbwork.stepsize();
          
          for_3 (i, j, k, bbwork, bbwork.stepsize())
            FastIndex3(tgt,i,j,k) = Processor();
          end_for
        } 
        
      int StepsizeMax = 1;
      for (int l=1; l<_Grid.NLevels(); l++)
        StepsizeMax *= _Grid.RefineFactor(l-1); 

      BBox bbox(DB().bbox());
      Coords step = DB().bbox().stepsize();
      Coords dimstep = DimUsed*step;
      bbox.upper() += dimstep;
      for_3 (i, j, k, bbox, step)
        v[idx+offset] = Processor() + 1.0;
      
        TheEvaluator().mvals[0] = Min(TheEvaluator().mvals[0], v[idx+offset]);
        TheEvaluator().mvals[1] = Max(TheEvaluator().mvals[1], v[idx+offset]);
        idx++;        
        end_for 

      EndFastIndex3(tgt); 
      offset+=idx;
    }

    else if (*key==2) {
      int idx=0;
      BBox lbbox(grow(DB().bbox(),DimUsed));
      data_block_type dbwork(lbbox);
      dbwork = -1.0;
      BeginFastIndex3(tgt, lbbox, dbwork.data(), DataType);
      
      for (gblock_list_iterator bit=_Grid.BlockList().begin();
           bit!=_Grid.BlockList().end(); bit++) 
        if ((*bit)->Level() == Level()) {
          BBox bbcur(growupper((*bit)->DB().bbox(),1));
          bbcur.setstepsize(lbbox.stepsize());
          bbcur.growupper(-1);
          BBox bbwork(lbbox * bbcur);
          if (bbwork.empty())
            continue;
          
          bbwork.lower() = (bbwork.lower()/bbwork.stepsize())*
                           bbwork.stepsize();
          bbwork.upper() = (bbwork.upper()/bbwork.stepsize())*
                           bbwork.stepsize();
          
          for_3 (i, j, k, bbwork, bbwork.stepsize())
            FastIndex3(tgt,i,j,k) = Level();
          end_for
        } 
        
      int StepsizeMax = 1;
      for (int l=1; l<_Grid.NLevels(); l++)
        StepsizeMax *= _Grid.RefineFactor(l-1); 

      BBox bbox(DB().bbox());
      Coords step = DB().bbox().stepsize();
      Coords dimstep = DimUsed*step;
      bbox.upper() += dimstep;
      for_3 (i, j, k, bbox, step)
        v[idx+offset] = Level() + 1.0;
      
        TheEvaluator().mvals[2] = Min(TheEvaluator().mvals[2], v[idx+offset]);
        TheEvaluator().mvals[3] = Max(TheEvaluator().mvals[3], v[idx+offset]);
        idx++;        
        end_for 

      EndFastIndex3(tgt); 
      offset+=idx;
    }

    else {
      BBox lbbox(grow(DB().bbox(),DimUsed));
      lbbox *= _Grid.GlobalBBox();
      lbbox.lower() = (lbbox.lower()/lbbox.stepsize())*lbbox.stepsize();
      lbbox.upper() = (lbbox.upper()/lbbox.stepsize())*lbbox.stepsize();
      vector_block_type dbwork(lbbox);
      dbwork = FLT_MAX;
      BeginFastIndex3(tgt, lbbox, dbwork.data(), VectorType);
        
      for (gblock_list_iterator bit=_Grid.BlockList().begin();
           bit!=_Grid.BlockList().end(); bit++) {
        BBox bbcur(growupper((*bit)->DB().bbox(),1));
        bbcur.setstepsize(lbbox.stepsize());
        bbcur.growupper(-1);
        BBox bbwork(lbbox * bbcur);
        if (bbwork.empty())
          continue;
        
        bbwork.lower() = (bbwork.lower()/bbwork.stepsize())*bbwork.stepsize();
        bbwork.upper() = (bbwork.upper()/bbwork.stepsize())*bbwork.stepsize();

        Coords step_src = (*bit)->DB().bbox().stepsize();
        Coords step_tgt = bbwork.stepsize();
        if (step_tgt(0) == step_src(0)) 
          dbwork.copy((*bit)->DB(), bbwork);
        
        else if (step_tgt(0) > step_src(0)) {
          BeginFastIndex3(src, (*bit)->DB().bbox(), (*bit)->DB().data(),
                          VectorType);
          for_3 (i, j, k, bbwork, bbwork.stepsize())
            BBox AddOver(3,i,j,k,
                         i+step_tgt(0)-step_src(0),
                         j+step_tgt(1)-step_src(1),
                         k+step_tgt(2)-step_src(2),
                         step_src(0),step_src(1),step_src(2));
            AddOver *= (*bit)->DB().bbox();
            float c=0;
            FastIndex3(tgt,i,j,k) = 0.0;
            for_3 (l, m, n, AddOver, step_src)
              FastIndex3(tgt,i,j,k) += FastIndex3(src,l,m,n);
              c += 1.0;
            end_for
            if (c>0) FastIndex3(tgt,i,j,k) /= c;
          end_for
          EndFastIndex3(src);
        }

        else {
          BeginFastIndex3(src, (*bit)->DB().bbox(), (*bit)->DB().data(),
                          VectorType);
          for_3 (i, j, k, bbwork, bbwork.stepsize())
            FastIndex3(tgt,i,j,k) = FastIndex3(src,i-i%step_src(0),
                                               j-j%step_src(1),
                                               k-k%step_src(2));
          end_for
          EndFastIndex3(src);
        }
      } 
      EndFastIndex3(tgt); 
        
      BBox bbox(DB().bbox());
      Coords dimstep = DimUsed*DB().bbox().stepsize();
      bbox.upper() += dimstep;
      TheEvaluator().SetScalNodes(*key,v,offset,dbwork,bbox,dx,CompRead);
    }
  }

  virtual void SetPointsCells(int *key,point_list_type &p,float v[],bool* CompRead) {
    register int d;
    float lb[3], ub[3];
    for (d=0; d<3; d++) {
      lb[d]=DB().lower(d)*dx[d]+geom[2*d];
      ub[d]=(DB().extents(d)+DimUsed(d)-1)*DB().stepsize(d)*dx[d]+lb[d];
    }
    std::vector<unsigned int> plist;
    for (register unsigned int i=0;i<p.size()/3;i++) {
      for (d=0; d<3; d++) 
        if (DimUsed(d))
          if (p[3*i+d]<lb[d] || p[3*i+d]>ub[d]) break;
      if (d==3) plist.push_back(i); 
    }

    if (!plist.empty()) {
      float *tmpdata = new float[NCells()];
      int offset=0;
      SetScalCells(key,tmpdata,offset,CompRead);
      for (register unsigned int i=0;i<plist.size();i++) {
        int idx[3];
        for (d=0; d<3; d++) {
          if (DimUsed(d))
            idx[d]=int((p[3*plist[i]+d]-lb[d])/(DB().stepsize(d)*dx[d]));
          else idx[d]=0;
          if (idx[d]<0 || idx[d]>=DB().extents(d)) break;
        }
        if (d==3) {
          v[plist[i]]=tmpdata[idx[0]+DB().extents(0)*(idx[1]+DB().extents(1)*idx[2])];
          
          // Do not test points anymore after data has been assigned once
          for (d=0; d<3; d++) 
            p[3*plist[i]+d]=-1.e37; 
        }
      }
      delete [] tmpdata;
    }
  }

  virtual void SetPointsNodes(int *key,point_list_type &p,float v[],bool* CompRead) {
    register int d;
    float lb[3], ub[3];
    for (d=0; d<3; d++) {
      lb[d]=DB().lower(d)*dx[d]+geom[2*d];
      ub[d]=(DB().extents(d)+DimUsed(d)-1)*DB().stepsize(d)*dx[d]+lb[d];
    }
    std::vector<unsigned int> plist;
    for (register unsigned int i=0;i<p.size()/3;i++) {
      for (d=0; d<3; d++) 
        if (DimUsed(d))
          if (p[3*i+d]<lb[d] || p[3*i+d]>ub[d]) break;
      if (d==3) plist.push_back(i); 
    }

    if (!plist.empty()) {
      float *tmpdata = new float[NNodes()];
      int offset=0;
      SetScalNodes(key,tmpdata,offset,CompRead);
      for (register unsigned int i=0;i<plist.size();i++) {
        int idx[3], idp[3], ext[3];
        float fac[3];
        for (d=0; d<3; d++) {
          if (DimUsed(d)) {
            idx[d]=int((p[3*plist[i]+d]-lb[d])/(DB().stepsize(d)*dx[d]));
            fac[d]=(p[3*plist[i]+d]-lb[d]-idx[d]*DB().stepsize(d)*dx[d])/
              (DB().stepsize(d)*dx[d]);
            if (fac[d]<=0.5) fac[d]=0.;
            else fac[d]=1.;
            idp[d]=idx[d]+1;
          }
          else {
            idx[d]=0; fac[d]=0.; idp[d]=0;
          }
          ext[d]=DB().extents(d)+DimUsed(d);
          if (idx[d]<0 || idp[d]>=ext[d]) break;
        }
        
        if (d==3) {
          float &d0=tmpdata[idx[0]+ext[0]*(idx[1]+ext[1]*idx[2])];
          float &d1=tmpdata[idp[0]+ext[0]*(idx[1]+ext[1]*idx[2])];
          float &d2=tmpdata[idx[0]+ext[0]*(idp[1]+ext[1]*idx[2])];
          float &d3=tmpdata[idp[0]+ext[0]*(idp[1]+ext[1]*idx[2])];
          float &d4=tmpdata[idx[0]+ext[0]*(idx[1]+ext[1]*idp[2])];
          float &d5=tmpdata[idp[0]+ext[0]*(idx[1]+ext[1]*idp[2])];
          float &d6=tmpdata[idx[0]+ext[0]*(idp[1]+ext[1]*idp[2])];
          float &d7=tmpdata[idp[0]+ext[0]*(idp[1]+ext[1]*idp[2])];
        
          float x0=FLT_MAX, x1=FLT_MAX, x2=FLT_MAX, x3=FLT_MAX;
          if (d0<FLT_MAX && d1<FLT_MAX) x0=(1.-fac[0])*d0+fac[0]*d1;
          else { if (d0<FLT_MAX) x0=d0; else if (d1<FLT_MAX) x0=d1; }
          if (d2<FLT_MAX && d3<FLT_MAX) x1=(1.-fac[0])*d2+fac[0]*d3;
          else { if (d2<FLT_MAX) x1=d2; else if (d3<FLT_MAX) x1=d3; }
          if (d4<FLT_MAX && d5<FLT_MAX) x2=(1.-fac[0])*d4+fac[0]*d5;
          else { if (d4<FLT_MAX) x2=d4; else if (d5<FLT_MAX) x2=d5; }
          if (d6<FLT_MAX && d7<FLT_MAX) x3=(1.-fac[0])*d6+fac[0]*d7;
          else { if (d6<FLT_MAX) x3=d6; else if (d7<FLT_MAX) x3=d7; }

          float y0=FLT_MAX, y1=FLT_MAX;
          if (x0<FLT_MAX && x1<FLT_MAX) y0=(1.-fac[1])*x0+fac[1]*x1;
          else { if (x0<FLT_MAX) y0=x0; else if (x1<FLT_MAX) y0=x1; }
          if (x2<FLT_MAX && x3<FLT_MAX) y1=(1.-fac[1])*x2+fac[1]*x3;
          else { if (x2<FLT_MAX) y1=x2; else if (x3<FLT_MAX) y1=x3; }

          v[plist[i]]=FLT_MAX;
          if (y0<FLT_MAX && y1<FLT_MAX) v[plist[i]]=(1.-fac[2])*y0+fac[2]*y1;
          else { if (y0<FLT_MAX) v[plist[i]]=y0; else if (y1<FLT_MAX) v[plist[i]]=y1; }

          // Do not test points anymore after data has been assigned once
          for (d=0; d<3; d++) 
            p[3*plist[i]+d]=-1.e37; 
        }
      }
      delete [] tmpdata;
    }
  }

  virtual void VectCells(int *key,float v[][3],int& offset,bool* CompRead) {
    std::cerr << ".";
    int idx=0;  
    BeginFastIndex3(dat, DB().bbox(), DB().data(), VectorType);
    for_3 (i,j,k,DB().bbox(),DB().bbox().stepsize())
      v[idx+offset][0]= FastIndex3(dat,i,j,k)(1);
      v[idx+offset][1]= FastIndex3(dat,i,j,k)(2);
      v[idx+offset][2]= FastIndex3(dat,i,j,k)(3);           
      idx++;          
    end_for
    offset+=idx;
    EndFastIndex3(dat); 
  }

  virtual void VectNodes(int *key,float v[][3],int& offset,bool* CompRead) {
    std::cerr << ".";
    int idx=0;

    BBox lbbox(grow(DB().bbox(),DimUsed));
    lbbox *= _Grid.GlobalBBox();
    lbbox.lower() = (lbbox.lower()/lbbox.stepsize())*lbbox.stepsize();
    lbbox.upper() = (lbbox.upper()/lbbox.stepsize())*lbbox.stepsize();
    vector_block_type dbwork(lbbox);
    dbwork = FLT_MAX;
    BeginFastIndex3(tgt, lbbox, dbwork.data(), VectorType);
    
    for (gblock_list_iterator bit=_Grid.BlockList().begin();
         bit!=_Grid.BlockList().end(); bit++) {
      BBox bbcur(growupper((*bit)->DB().bbox(),1));
      bbcur.setstepsize(lbbox.stepsize());
      bbcur.growupper(-1);
      BBox bbwork(lbbox * bbcur);
      if (bbwork.empty())
        continue;

      bbwork.lower() = (bbwork.lower()/bbwork.stepsize())*bbwork.stepsize();
      bbwork.upper() = (bbwork.upper()/bbwork.stepsize())*bbwork.stepsize();
      
      Coords step_src = (*bit)->DB().bbox().stepsize();
      Coords step_tgt = bbwork.stepsize();
      if (step_tgt(0) == step_src(0)) 
        dbwork.copy((*bit)->DB(), bbwork);

      else if (step_tgt(0) > step_src(0)) {
        BeginFastIndex3(src, (*bit)->DB().bbox(), (*bit)->DB().data(),
                        VectorType);
        for_3 (i, j, k, bbwork, bbwork.stepsize())
          BBox AddOver(3,i,j,k,
                       i+step_tgt(0)-step_src(0),
                       j+step_tgt(1)-step_src(1),
                       k+step_tgt(2)-step_src(2),
                       step_src(0),step_src(1),step_src(2));
          AddOver *= (*bit)->DB().bbox();
          float c=0;
          FastIndex3(tgt,i,j,k) = 0.0;
          for_3 (l, m, n, AddOver, step_src)
            FastIndex3(tgt,i,j,k) += FastIndex3(src,l,m,n);
            c += 1.0;
          end_for
          if (c>0) FastIndex3(tgt,i,j,k) /= c;
        end_for
        EndFastIndex3(src);
      }

      else {
        BeginFastIndex3(src, (*bit)->DB().bbox(), (*bit)->DB().data(),
                        VectorType);
        for_3 (i, j, k, bbwork, bbwork.stepsize())
          FastIndex3(tgt,i,j,k) = FastIndex3(src,i-i%step_src(0),
                                             j-j%step_src(1),
                                             k-k%step_src(2));
        end_for
        EndFastIndex3(src);
      }
    } 
      
    BBox bbox(DB().bbox());
    Coords dimstep = DimUsed*DB().bbox().stepsize();
    bbox.upper() += dimstep;
      
    for_3 (i, j, k, bbox, dimstep)
      BBox AddOver(3,i-dimstep(0),j-dimstep(1),k-dimstep(2),
                   i,j,k,dimstep(0),dimstep(1),dimstep(2));
      AddOver *= lbbox;
          
      float c=0;
      v[idx+offset][0]=0;
      v[idx+offset][1]=0;
      v[idx+offset][2]=0;
      for_3 (l, m, n, AddOver, dimstep)
        v[idx+offset][0] += FastIndex3(tgt,l,m,n)(1);
        v[idx+offset][1] += FastIndex3(tgt,l,m,n)(2);
        v[idx+offset][2] += FastIndex3(tgt,l,m,n)(3);
        c += 1.0;
      end_for
      if (c>0) {
        v[idx+offset][0] /= c;
        v[idx+offset][1] /= c;
        v[idx+offset][2] /= c;
      }
      idx++;          
    end_for 
      
    EndFastIndex3(tgt); 
    offset+=idx;
  }

  virtual vector_block_type& DB() { return *_data_block; }
  virtual int NNodes() { return (!DB().bbox().empty() ? 
     ((DB().extents(0)+DimUsed(0))*
      (DB().extents(1)+DimUsed(1))*
      (DB().extents(2)+DimUsed(2))) : 0); }
  virtual int NCells() { return DB().size(); }
  virtual const int& Level() const { return _Level; } 
  virtual int Level() { return _Level; } 
  virtual const int& Processor() const { return _Processor; } 
  virtual int Processor() { return _Processor; } 
  virtual evaluator_type& TheEvaluator() { return _Evaluator; }

protected:
  grid_type& _Grid;
  vector_block_type* _data_block;
  evaluator_type& _Evaluator;
  int _Level;
  int _Processor;
  float* dx;
  float* geom;
  int _nodeOffset;
  Coords DimUsed;
};

#endif

---