00001
00002
00003
00004
00005
00006
00007
00008
00009 #ifndef AMROC_VISUALGRIDDATA_H
00010 #define AMROC_VISUALGRIDDATA_H
00011
00019 #include "GridData3.h"
00020 #include "Evaluator.h"
00021
00022 #include <iostream>
00023 #include <vector>
00024 #include <cfloat>
00025 #include <cmath>
00026
00027 template <class VectorType> class VisualGrid;
00028
00035 template <class VectorType>
00036 class VisualGridData {
00037 typedef typename VectorType::InternalDataType DataType;
00038 typedef VisualGridData<VectorType> self;
00039
00040 public:
00041 typedef VisualGrid<VectorType> grid_type;
00042 typedef Evaluator<VectorType> evaluator_type;
00043 typedef GridData<DataType,3> data_block_type;
00044 typedef GridData<VectorType,3> vector_block_type;
00045
00046 typedef std::vector<self *> gblock_list_type;
00047 typedef typename gblock_list_type::iterator gblock_list_iterator;
00048 typedef std::vector<float> point_list_type;
00049 typedef typename point_list_type::iterator point_list_iterator;
00050
00051 enum direction {x_dir, y_dir, z_dir};
00052
00053 public:
00054 VisualGridData(grid_type& gr, data_block_type& db, BBox &bb, int comp, evaluator_type& ev,
00055 int lev, int proc, float* pdx, float* pgeom) :
00056 _Grid(gr), _Evaluator(ev), _Level(lev), _Processor(proc), dx(pdx), geom(pgeom) {
00057
00058 _data_block = new vector_block_type(bb);
00059 BeginFastIndex3(source, db.bbox(), db.data(), DataType);
00060 BeginFastIndex3(target, _data_block->bbox(), _data_block->data(),
00061 VectorType);
00062 for_3 (i, j, k, bb, bb.stepsize())
00063 FastIndex3(target,i,j,k)(comp) = FastIndex3(source,i,j,k);
00064 end_for
00065 EndFastIndex3(source);
00066 EndFastIndex3(target);
00067 DimUsed = gr.Dimensions();
00068 }
00069
00070 VisualGridData(grid_type& gr, vector_block_type& db, BBox &bb, evaluator_type& ev,
00071 int lev, int proc, float* pdx, float* pgeom) :
00072 _Grid(gr), _Evaluator(ev), _Level(lev), _Processor(proc), dx(pdx), geom(pgeom) {
00073
00074 _data_block = new vector_block_type(bb);
00075 BeginFastIndex3(source, db.bbox(), db.data(), VectorType);
00076 BeginFastIndex3(target, _data_block->bbox(), _data_block->data(),
00077 VectorType);
00078 for_3 (i, j, k, bb, bb.stepsize())
00079 FastIndex3(target,i,j,k) = FastIndex3(source,i,j,k);
00080 end_for
00081 EndFastIndex3(source);
00082 EndFastIndex3(target);
00083 DimUsed = gr.Dimensions();
00084 }
00085
00086 virtual ~VisualGridData()
00087 { delete _data_block; }
00088
00089 public:
00090 virtual void SetBlocks(int blocks[],int& offset,int& nOffset) {
00091 std::cerr << ".";
00092 _nodeOffset = nOffset;
00093 blocks[offset] = DB().extents(0)+DimUsed(0);
00094 blocks[offset+1] = DB().extents(1)+DimUsed(1);
00095 blocks[offset+2] = DB().extents(2)+DimUsed(2);
00096 nOffset += blocks[offset]*blocks[offset+1]*blocks[offset+2];
00097 offset += 3;
00098 }
00099
00100 virtual void SetGridCells(float xyz[][3],int& offset){
00101 std::cerr << ".";
00102 int i, j, k,orii,orij,orik,ssi,ssj,ssk,ijk,idx;
00103 int exti,extj,extk;
00104 exti=DB().extents(0);
00105 extj=DB().extents(1);
00106 extk=DB().extents(2);
00107 orii=DB().lower(0);
00108 orij=DB().lower(1);
00109 orik=DB().lower(2);
00110 ssi=DB().stepsize(0);
00111 ssj=DB().stepsize(1);
00112 ssk=DB().stepsize(2);
00113
00114 idx = 0;
00115 for (k = 0; k < extk; k++){
00116 for (j = 0; j < extj; j++){
00117 for (i = 0; i < exti; i++){
00118
00119 ijk = k*extj*exti + j*exti + i;
00120
00121 xyz[ijk+offset][0] = (orii+((float)i+0.5)*ssi)*dx[0]+geom[0];
00122 xyz[ijk+offset][1] = (orij+((float)j+0.5)*ssj)*dx[1]+geom[2];
00123 xyz[ijk+offset][2] = (orik+((float)k+0.5)*ssk)*dx[2]+geom[4];
00124 idx++;
00125 }
00126 }
00127 }
00128 offset+=idx;
00129 }
00130 virtual void SetGridNodes(float xyz[][3],int& offset) {
00131 std::cerr << ".";
00132 int i, j, k,orii,orij,orik,ssi,ssj,ssk,ijk,idx;
00133 int exti,extj,extk;
00134 exti=DB().extents(0)+DimUsed(0);
00135 extj=DB().extents(1)+DimUsed(1);
00136 extk=DB().extents(2)+DimUsed(2);
00137 orii=DB().lower(0);
00138 orij=DB().lower(1);
00139 orik=DB().lower(2);
00140 ssi=DB().stepsize(0);
00141 ssj=DB().stepsize(1);
00142 ssk=DB().stepsize(2);
00143
00144 idx = 0;
00145 for (k = 0; k < extk; k++){
00146 for (j = 0; j < extj; j++){
00147 for (i = 0; i < exti; i++){
00148
00149 ijk = k*extj*exti + j*exti + i;
00150
00151 xyz[ijk+offset][0] = (float)(orii+i*ssi)*dx[0]+geom[0];
00152 xyz[ijk+offset][1] = (float)(orij+j*ssj)*dx[1]+geom[2];
00153 xyz[ijk+offset][2] = (float)(orik+k*ssk)*dx[2]+geom[4];
00154 idx++;
00155 }
00156 }
00157 }
00158 offset+=idx;
00159 }
00160
00161 virtual void SetGridConns(int con[][8],int& offset) {
00162 std::cerr << ".";
00163 int i,j,k,ijk,idx;
00164 int exti,extj,extk;
00165 exti=DB().extents(0)+DimUsed(0);
00166 extj=DB().extents(1)+DimUsed(1);
00167 extk=DB().extents(2)+DimUsed(2);
00168
00169 idx = 0;
00170 for (k = 0; k < extk-DimUsed(2); k++){
00171 for (j = 0; j < extj-DimUsed(1); j++){
00172 for (i = 0; i < exti-DimUsed(0); i++){
00173
00174 ijk = k*(extj-DimUsed(1))*(exti-DimUsed(0)) +
00175 j*(exti-DimUsed(0)) + i;
00176
00177 con[ijk+offset][0] = k*extj*exti +
00178 j*exti +
00179 i + _nodeOffset;
00180 con[ijk+offset][1] = k*extj*exti +
00181 j*exti +
00182 (i+DimUsed(0)) + _nodeOffset;
00183 con[ijk+offset][2] = k*extj*exti +
00184 (j+DimUsed(1))*exti +
00185 i + _nodeOffset;
00186 con[ijk+offset][3] = k*extj*exti +
00187 (j+DimUsed(1))*exti +
00188 (i+DimUsed(0)) + _nodeOffset;
00189 con[ijk+offset][4] = (k+DimUsed(2))*extj*exti +
00190 j*exti +
00191 i + _nodeOffset;
00192 con[ijk+offset][5] = (k+DimUsed(2))*extj*exti +
00193 j*exti +
00194 (i+DimUsed(0)) + _nodeOffset;
00195 con[ijk+offset][6] = (k+DimUsed(2))*extj*exti +
00196 (j+DimUsed(1))*exti +
00197 i + _nodeOffset;
00198 con[ijk+offset][7] = (k+DimUsed(2))*extj*exti +
00199 (j+DimUsed(1))*exti +
00200 (i+DimUsed(0)) + _nodeOffset;
00201
00202 idx++;
00203 }
00204 }
00205 }
00206 offset+=idx;
00207 }
00208
00209 virtual void SetScalCells(int *key,float v[],int& offset,bool* CompRead) {
00210 std::cerr << ".";
00211 if (*key==1 || *key==2) {
00212 int idx=0;
00213 for_3 (i,j,k,DB().bbox(),DB().bbox().stepsize())
00214 if (*key==1) {
00215 v[idx+offset]= Processor() + 1.0;
00216 TheEvaluator().mvals[0] = Min(TheEvaluator().mvals[0],v[idx+offset]);
00217 TheEvaluator().mvals[1] = Max(TheEvaluator().mvals[1],v[idx+offset]);
00218 }
00219 else {
00220 v[idx+offset]= Level() + 1.0;
00221 TheEvaluator().mvals[2] = Min(TheEvaluator().mvals[2],v[idx+offset]);
00222 TheEvaluator().mvals[3] = Max(TheEvaluator().mvals[3],v[idx+offset]);
00223 }
00224 idx++;
00225 end_for
00226 offset+=idx;
00227 }
00228 else
00229 TheEvaluator().SetScalCells(*key,v,offset,DB(),dx,CompRead);
00230 }
00231
00232 virtual void SetScalNodes(int *key,float v[],int& offset,bool* CompRead) {
00233 std::cerr << ".";
00234 if (*key==1) {
00235 int idx=0;
00236 BBox lbbox(grow(DB().bbox(),DimUsed));
00237 data_block_type dbwork(lbbox);
00238 dbwork = -1.0;
00239 BeginFastIndex3(tgt, lbbox, dbwork.data(), DataType);
00240
00241 for (gblock_list_iterator bit=_Grid.BlockList().begin();
00242 bit!=_Grid.BlockList().end(); bit++)
00243 if ((*bit)->Processor() == Processor()) {
00244 BBox bbcur(growupper((*bit)->DB().bbox(),1));
00245 bbcur.setstepsize(lbbox.stepsize());
00246 bbcur.growupper(-1);
00247 BBox bbwork(lbbox * bbcur);
00248 if (bbwork.empty())
00249 continue;
00250
00251 bbwork.lower() = (bbwork.lower()/bbwork.stepsize())*
00252 bbwork.stepsize();
00253 bbwork.upper() = (bbwork.upper()/bbwork.stepsize())*
00254 bbwork.stepsize();
00255
00256 for_3 (i, j, k, bbwork, bbwork.stepsize())
00257 FastIndex3(tgt,i,j,k) = Processor();
00258 end_for
00259 }
00260
00261 int StepsizeMax = 1;
00262 for (int l=1; l<_Grid.NLevels(); l++)
00263 StepsizeMax *= _Grid.RefineFactor(l-1);
00264
00265 BBox bbox(DB().bbox());
00266 Coords step = DB().bbox().stepsize();
00267 Coords dimstep = DimUsed*step;
00268 bbox.upper() += dimstep;
00269 for_3 (i, j, k, bbox, step)
00270 v[idx+offset] = Processor() + 1.0;
00271
00272 TheEvaluator().mvals[0] = Min(TheEvaluator().mvals[0], v[idx+offset]);
00273 TheEvaluator().mvals[1] = Max(TheEvaluator().mvals[1], v[idx+offset]);
00274 idx++;
00275 end_for
00276
00277 EndFastIndex3(tgt);
00278 offset+=idx;
00279 }
00280
00281 else if (*key==2) {
00282 int idx=0;
00283 BBox lbbox(grow(DB().bbox(),DimUsed));
00284 data_block_type dbwork(lbbox);
00285 dbwork = -1.0;
00286 BeginFastIndex3(tgt, lbbox, dbwork.data(), DataType);
00287
00288 for (gblock_list_iterator bit=_Grid.BlockList().begin();
00289 bit!=_Grid.BlockList().end(); bit++)
00290 if ((*bit)->Level() == Level()) {
00291 BBox bbcur(growupper((*bit)->DB().bbox(),1));
00292 bbcur.setstepsize(lbbox.stepsize());
00293 bbcur.growupper(-1);
00294 BBox bbwork(lbbox * bbcur);
00295 if (bbwork.empty())
00296 continue;
00297
00298 bbwork.lower() = (bbwork.lower()/bbwork.stepsize())*
00299 bbwork.stepsize();
00300 bbwork.upper() = (bbwork.upper()/bbwork.stepsize())*
00301 bbwork.stepsize();
00302
00303 for_3 (i, j, k, bbwork, bbwork.stepsize())
00304 FastIndex3(tgt,i,j,k) = Level();
00305 end_for
00306 }
00307
00308 int StepsizeMax = 1;
00309 for (int l=1; l<_Grid.NLevels(); l++)
00310 StepsizeMax *= _Grid.RefineFactor(l-1);
00311
00312 BBox bbox(DB().bbox());
00313 Coords step = DB().bbox().stepsize();
00314 Coords dimstep = DimUsed*step;
00315 bbox.upper() += dimstep;
00316 for_3 (i, j, k, bbox, step)
00317 v[idx+offset] = Level() + 1.0;
00318
00319 TheEvaluator().mvals[2] = Min(TheEvaluator().mvals[2], v[idx+offset]);
00320 TheEvaluator().mvals[3] = Max(TheEvaluator().mvals[3], v[idx+offset]);
00321 idx++;
00322 end_for
00323
00324 EndFastIndex3(tgt);
00325 offset+=idx;
00326 }
00327
00328 else {
00329 BBox lbbox(grow(DB().bbox(),DimUsed));
00330 lbbox *= _Grid.GlobalBBox();
00331 lbbox.lower() = (lbbox.lower()/lbbox.stepsize())*lbbox.stepsize();
00332 lbbox.upper() = (lbbox.upper()/lbbox.stepsize())*lbbox.stepsize();
00333 vector_block_type dbwork(lbbox);
00334 dbwork = FLT_MAX;
00335 BeginFastIndex3(tgt, lbbox, dbwork.data(), VectorType);
00336
00337 for (gblock_list_iterator bit=_Grid.BlockList().begin();
00338 bit!=_Grid.BlockList().end(); bit++) {
00339 BBox bbcur(growupper((*bit)->DB().bbox(),1));
00340 bbcur.setstepsize(lbbox.stepsize());
00341 bbcur.growupper(-1);
00342 BBox bbwork(lbbox * bbcur);
00343 if (bbwork.empty())
00344 continue;
00345
00346 bbwork.lower() = (bbwork.lower()/bbwork.stepsize())*bbwork.stepsize();
00347 bbwork.upper() = (bbwork.upper()/bbwork.stepsize())*bbwork.stepsize();
00348
00349 Coords step_src = (*bit)->DB().bbox().stepsize();
00350 Coords step_tgt = bbwork.stepsize();
00351 if (step_tgt(0) == step_src(0))
00352 dbwork.copy((*bit)->DB(), bbwork);
00353
00354 else if (step_tgt(0) > step_src(0)) {
00355 BeginFastIndex3(src, (*bit)->DB().bbox(), (*bit)->DB().data(),
00356 VectorType);
00357 for_3 (i, j, k, bbwork, bbwork.stepsize())
00358 BBox AddOver(3,i,j,k,
00359 i+step_tgt(0)-step_src(0),
00360 j+step_tgt(1)-step_src(1),
00361 k+step_tgt(2)-step_src(2),
00362 step_src(0),step_src(1),step_src(2));
00363 AddOver *= (*bit)->DB().bbox();
00364 float c=0;
00365 FastIndex3(tgt,i,j,k) = 0.0;
00366 for_3 (l, m, n, AddOver, step_src)
00367 FastIndex3(tgt,i,j,k) += FastIndex3(src,l,m,n);
00368 c += 1.0;
00369 end_for
00370 if (c>0) FastIndex3(tgt,i,j,k) /= c;
00371 end_for
00372 EndFastIndex3(src);
00373 }
00374
00375 else {
00376 BeginFastIndex3(src, (*bit)->DB().bbox(), (*bit)->DB().data(),
00377 VectorType);
00378 for_3 (i, j, k, bbwork, bbwork.stepsize())
00379 FastIndex3(tgt,i,j,k) = FastIndex3(src,i-i%step_src(0),
00380 j-j%step_src(1),
00381 k-k%step_src(2));
00382 end_for
00383 EndFastIndex3(src);
00384 }
00385 }
00386 EndFastIndex3(tgt);
00387
00388 BBox bbox(DB().bbox());
00389 Coords dimstep = DimUsed*DB().bbox().stepsize();
00390 bbox.upper() += dimstep;
00391 TheEvaluator().SetScalNodes(*key,v,offset,dbwork,bbox,dx,CompRead);
00392 }
00393 }
00394
00395 virtual void SetPointsCells(int *key,point_list_type &p,float v[],bool* CompRead) {
00396 register int d;
00397 float lb[3], ub[3];
00398 for (d=0; d<3; d++) {
00399 lb[d]=DB().lower(d)*dx[d]+geom[2*d];
00400 ub[d]=(DB().extents(d)+DimUsed(d)-1)*DB().stepsize(d)*dx[d]+lb[d];
00401 }
00402 std::vector<unsigned int> plist;
00403 for (register unsigned int i=0;i<p.size()/3;i++) {
00404 for (d=0; d<3; d++)
00405 if (DimUsed(d))
00406 if (p[3*i+d]<lb[d] || p[3*i+d]>ub[d]) break;
00407 if (d==3) plist.push_back(i);
00408 }
00409
00410 if (!plist.empty()) {
00411 float *tmpdata = new float[NCells()];
00412 int offset=0;
00413 SetScalCells(key,tmpdata,offset,CompRead);
00414 for (register unsigned int i=0;i<plist.size();i++) {
00415 int idx[3];
00416 for (d=0; d<3; d++) {
00417 if (DimUsed(d))
00418 idx[d]=int((p[3*plist[i]+d]-lb[d])/(DB().stepsize(d)*dx[d]));
00419 else idx[d]=0;
00420 assert (idx[d]>=0 && idx[d]<DB().extents(d));
00421
00422 p[3*plist[i]+d]=-1.e37;
00423 }
00424 v[plist[i]]=tmpdata[idx[0]+DB().extents(0)*(idx[1]+DB().extents(1)*idx[2])];
00425 }
00426 delete [] tmpdata;
00427 }
00428 }
00429
00430 virtual void SetPointsNodes(int *key,point_list_type &p,float v[],bool* CompRead) {
00431 register int d;
00432 float lb[3], ub[3];
00433 for (d=0; d<3; d++) {
00434 lb[d]=DB().lower(d)*dx[d]+geom[2*d];
00435 ub[d]=(DB().extents(d)+DimUsed(d)-1)*DB().stepsize(d)*dx[d]+lb[d];
00436 }
00437 std::vector<unsigned int> plist;
00438 for (register unsigned int i=0;i<p.size()/3;i++) {
00439 for (d=0; d<3; d++)
00440 if (DimUsed(d))
00441 if (p[3*i+d]<lb[d] || p[3*i+d]>ub[d]) break;
00442 if (d==3) plist.push_back(i);
00443 }
00444
00445 if (!plist.empty()) {
00446 float *tmpdata = new float[NNodes()];
00447 int offset=0;
00448 SetScalNodes(key,tmpdata,offset,CompRead);
00449 for (register unsigned int i=0;i<plist.size();i++) {
00450 int idx[3], idp[3], ext[3];
00451 float fac[3];
00452 for (d=0; d<3; d++) {
00453 if (DimUsed(d)) {
00454 idx[d]=int((p[3*plist[i]+d]-lb[d])/(DB().stepsize(d)*dx[d]));
00455 fac[d]=(p[3*plist[i]+d]-lb[d]-idx[d]*DB().stepsize(d)*dx[d])/
00456 (DB().stepsize(d)*dx[d]);
00457 idp[d]=idx[d]+1;
00458 }
00459 else {
00460 idx[d]=0; fac[d]=0.; idp[d]=0;
00461 }
00462 ext[d]=DB().extents(d)+DimUsed(d);
00463 assert (idx[d]>=0 && idp[d]<ext[d]);
00464
00465 p[3*plist[i]+d]=-1.e37;
00466 }
00467
00468 float &d0=tmpdata[idx[0]+ext[0]*(idx[1]+ext[1]*idx[2])];
00469 float &d1=tmpdata[idp[0]+ext[0]*(idx[1]+ext[1]*idx[2])];
00470 float &d2=tmpdata[idx[0]+ext[0]*(idp[1]+ext[1]*idx[2])];
00471 float &d3=tmpdata[idp[0]+ext[0]*(idp[1]+ext[1]*idx[2])];
00472 float &d4=tmpdata[idx[0]+ext[0]*(idx[1]+ext[1]*idp[2])];
00473 float &d5=tmpdata[idp[0]+ext[0]*(idx[1]+ext[1]*idp[2])];
00474 float &d6=tmpdata[idx[0]+ext[0]*(idp[1]+ext[1]*idp[2])];
00475 float &d7=tmpdata[idp[0]+ext[0]*(idp[1]+ext[1]*idp[2])];
00476
00477 v[plist[i]]=(1.0-fac[2])*((1.-fac[1])*((1.-fac[0])*d0+fac[0]*d1)+
00478 fac[1] *((1.-fac[0])*d2+fac[0]*d3))+
00479 fac[2] *((1.-fac[1])*((1.-fac[0])*d4+fac[0]*d5)+
00480 fac[1] *((1.-fac[0])*d6+fac[0]*d7));
00481 }
00482 delete [] tmpdata;
00483 }
00484 }
00485
00486 virtual void VectCells(int *key,float v[][3],int& offset,bool* CompRead) {
00487 std::cerr << ".";
00488 int idx=0;
00489 BeginFastIndex3(dat, DB().bbox(), DB().data(), VectorType);
00490 for_3 (i,j,k,DB().bbox(),DB().bbox().stepsize())
00491 v[idx+offset][0]= FastIndex3(dat,i,j,k)(1);
00492 v[idx+offset][1]= FastIndex3(dat,i,j,k)(2);
00493 v[idx+offset][2]= FastIndex3(dat,i,j,k)(3);
00494 idx++;
00495 end_for
00496 offset+=idx;
00497 EndFastIndex3(dat);
00498 }
00499
00500 virtual void VectNodes(int *key,float v[][3],int& offset,bool* CompRead) {
00501 std::cerr << ".";
00502 int idx=0;
00503
00504 BBox lbbox(grow(DB().bbox(),DimUsed));
00505 lbbox *= _Grid.GlobalBBox();
00506 lbbox.lower() = (lbbox.lower()/lbbox.stepsize())*lbbox.stepsize();
00507 lbbox.upper() = (lbbox.upper()/lbbox.stepsize())*lbbox.stepsize();
00508 vector_block_type dbwork(lbbox);
00509 dbwork = FLT_MAX;
00510 BeginFastIndex3(tgt, lbbox, dbwork.data(), VectorType);
00511
00512 for (gblock_list_iterator bit=_Grid.BlockList().begin();
00513 bit!=_Grid.BlockList().end(); bit++) {
00514 BBox bbcur(growupper((*bit)->DB().bbox(),1));
00515 bbcur.setstepsize(lbbox.stepsize());
00516 bbcur.growupper(-1);
00517 BBox bbwork(lbbox * bbcur);
00518 if (bbwork.empty())
00519 continue;
00520
00521 bbwork.lower() = (bbwork.lower()/bbwork.stepsize())*bbwork.stepsize();
00522 bbwork.upper() = (bbwork.upper()/bbwork.stepsize())*bbwork.stepsize();
00523
00524 Coords step_src = (*bit)->DB().bbox().stepsize();
00525 Coords step_tgt = bbwork.stepsize();
00526 if (step_tgt(0) == step_src(0))
00527 dbwork.copy((*bit)->DB(), bbwork);
00528
00529 else if (step_tgt(0) > step_src(0)) {
00530 BeginFastIndex3(src, (*bit)->DB().bbox(), (*bit)->DB().data(),
00531 VectorType);
00532 for_3 (i, j, k, bbwork, bbwork.stepsize())
00533 BBox AddOver(3,i,j,k,
00534 i+step_tgt(0)-step_src(0),
00535 j+step_tgt(1)-step_src(1),
00536 k+step_tgt(2)-step_src(2),
00537 step_src(0),step_src(1),step_src(2));
00538 AddOver *= (*bit)->DB().bbox();
00539 float c=0;
00540 FastIndex3(tgt,i,j,k) = 0.0;
00541 for_3 (l, m, n, AddOver, step_src)
00542 FastIndex3(tgt,i,j,k) += FastIndex3(src,l,m,n);
00543 c += 1.0;
00544 end_for
00545 if (c>0) FastIndex3(tgt,i,j,k) /= c;
00546 end_for
00547 EndFastIndex3(src);
00548 }
00549
00550 else {
00551 BeginFastIndex3(src, (*bit)->DB().bbox(), (*bit)->DB().data(),
00552 VectorType);
00553 for_3 (i, j, k, bbwork, bbwork.stepsize())
00554 FastIndex3(tgt,i,j,k) = FastIndex3(src,i-i%step_src(0),
00555 j-j%step_src(1),
00556 k-k%step_src(2));
00557 end_for
00558 EndFastIndex3(src);
00559 }
00560 }
00561
00562 BBox bbox(DB().bbox());
00563 Coords dimstep = DimUsed*DB().bbox().stepsize();
00564 bbox.upper() += dimstep;
00565
00566 for_3 (i, j, k, bbox, dimstep)
00567 BBox AddOver(3,i-dimstep(0),j-dimstep(1),k-dimstep(2),
00568 i,j,k,dimstep(0),dimstep(1),dimstep(2));
00569 AddOver *= lbbox;
00570
00571 float c=0;
00572 v[idx+offset][0]=0;
00573 v[idx+offset][1]=0;
00574 v[idx+offset][2]=0;
00575 for_3 (l, m, n, AddOver, dimstep)
00576 v[idx+offset][0] += FastIndex3(tgt,l,m,n)(1);
00577 v[idx+offset][1] += FastIndex3(tgt,l,m,n)(2);
00578 v[idx+offset][2] += FastIndex3(tgt,l,m,n)(3);
00579 c += 1.0;
00580 end_for
00581 if (c>0) {
00582 v[idx+offset][0] /= c;
00583 v[idx+offset][1] /= c;
00584 v[idx+offset][2] /= c;
00585 }
00586 idx++;
00587 end_for
00588
00589 EndFastIndex3(tgt);
00590 offset+=idx;
00591 }
00592
00593 virtual vector_block_type& DB() { return *_data_block; }
00594 virtual int NNodes() { return (!DB().bbox().empty() ?
00595 ((DB().extents(0)+DimUsed(0))*
00596 (DB().extents(1)+DimUsed(1))*
00597 (DB().extents(2)+DimUsed(2))) : 0); }
00598 virtual int NCells() { return DB().size(); }
00599 virtual const int& Level() const { return _Level; }
00600 virtual int Level() { return _Level; }
00601 virtual const int& Processor() const { return _Processor; }
00602 virtual int Processor() { return _Processor; }
00603 virtual evaluator_type& TheEvaluator() { return _Evaluator; }
00604
00605 protected:
00606 grid_type& _Grid;
00607 vector_block_type* _data_block;
00608 evaluator_type& _Evaluator;
00609 int _Level;
00610 int _Processor;
00611 float* dx;
00612 float* geom;
00613 int _nodeOffset;
00614 Coords DimUsed;
00615 };
00616
00617 #endif