00001
00002
00003 #ifndef AMROC_V3VISUALGRIDDATA_H
00004 #define AMROC_V3VISUALGRIDDATA_H
00005
00013 #include "VisualGridData.h"
00014
00015 float tol=FLT_EPSILON*100.0;
00016
00017 template <class VectorType> class V3VisualGrid;
00018
00025 template <class VectorType>
00026 class V3VisualGridData : public VisualGridData<VectorType> {
00027 typedef VisualGridData<VectorType> base;
00028 typedef V3VisualGridData<VectorType> self;
00029 typedef typename VectorType::InternalDataType DataType;
00030
00031 public:
00032 typedef V3VisualGrid<VectorType> grid_type;
00033 typedef typename base::grid_type base_grid_type;
00034 typedef typename base::evaluator_type evaluator_type;
00035 typedef typename base::data_block_type data_block_type;
00036 typedef typename base::vector_block_type vector_block_type;
00037
00038 V3VisualGridData(grid_type& gr, data_block_type& db, BBox &bb, int comp, evaluator_type& ev,
00039 int lev, int proc, float* pdx, float* pgeom) :
00040 base((base_grid_type &) gr, db, bb, comp, ev, lev, proc, pdx, pgeom) {
00041 CrossPlane = (float *)0;
00042 }
00043
00044 V3VisualGridData(grid_type& gr, vector_block_type& db, BBox &bb, evaluator_type& ev,
00045 int lev, int proc, float* pdx, float* pgeom) :
00046 base((base_grid_type &) gr, db, bb, ev, lev, proc, pdx, pgeom) {
00047 CrossPlane = (float *)0;
00048 }
00049
00050 void Draw3DGrids()
00051 {
00052 float orii,orij,orik;
00053 float exti,extj,extk;
00054 exti=base::DB().stepsize(0)*base::DB().extents(0)*base::dx[0];
00055 extj=base::DB().stepsize(1)*base::DB().extents(1)*base::dx[1];
00056 extk=base::DB().stepsize(2)*base::DB().extents(2)*base::dx[2];
00057
00058 orii=base::DB().lower(0)*base::dx[0]+base::geom[0];
00059 orij=base::DB().lower(1)*base::dx[1]+base::geom[2];
00060 orik=base::DB().lower(2)*base::dx[2]+base::geom[4];
00061
00062 float radii = 0.;
00063 int intzero = 0;
00064 int itype = 4;
00065 int icolor = 1;
00066 int np = 5;
00067 float color;
00068 switch (base::_Level){
00069 case 0:
00070 color=.54;
00071 break;
00072 case 1:
00073 color=.44;
00074 break;
00075 case 2:
00076 color=.34;
00077 break;
00078 case 3:
00079 color=.24;
00080 break;
00081 case 4:
00082 color=.14;
00083 break;
00084 default:
00085 color=.04;
00086 break;
00087 }
00088
00089
00090 long iv[1]={1};
00091 float rv[3]={color,color,color};
00092 float saveentry[3];
00093 int opt=18;
00094 char blank=' ';
00095 V3_GETSTATE(&opt,iv,saveentry,&blank,1);
00096 V3_SETSTATE(&opt,iv,rv,&blank,1);
00097
00098 itype=4;
00099 np=5;
00100 float col[8]={0.,0.,0.,0.,0.,0.,0.,0.};
00101 float line1[5*3] = { orii, orij, orik, orii, orij+extj, orik,
00102 orii+exti, orij+extj, orik, orii+exti, orij, orik,
00103 orii, orij, orik};
00104
00105 V3_OBJECT3D(&itype, &icolor, line1, &radii, col, &intzero, &np);
00106
00107 float line2[5*3] = { orii, orij, orik+extk,
00108 orii, orij+extj, orik+extk,
00109 orii+exti, orij+extj, orik+extk,
00110 orii+exti, orij, orik+extk,
00111 orii, orij, orik+extk};
00112
00113 V3_OBJECT3D(&itype, &icolor, line2, &radii, col, &intzero, &np);
00114
00115 itype = 3;
00116 np = 2*4;
00117 float lines[2*4*3] = { orii, orij, orik,
00118 orii, orij, orik+extk,
00119 orii, orij+extj, orik,
00120 orii, orij+extj, orik+extk,
00121
00122 orii+exti, orij+extj, orik,
00123 orii+exti, orij+extj, orik+extk,
00124 orii+exti, orij, orik,
00125 orii+exti, orij, orik+extk};
00126
00127 V3_OBJECT3D(&itype, &icolor, lines, &radii, col, &intzero, &np);
00128 V3_SETSTATE(&opt,iv,saveentry,&blank,1);
00129 }
00130
00131 void Draw3DCrossSections(){
00132 if (CrossPlane == (float *)0)
00133 return;
00134
00135 float ext[3]={base::DB().stepsize(0)*base::DB().extents(0)*base::dx[0],
00136 base::DB().stepsize(1)*base::DB().extents(1)*base::dx[1],
00137 base::DB().stepsize(2)*base::DB().extents(2)*base::dx[2]};
00138 float ori[3]={base::DB().lower(0)*base::dx[0]+base::geom[0],
00139 base::DB().lower(1)*base::dx[1]+base::geom[2],
00140 base::DB().lower(2)*base::dx[2]+base::geom[4]};
00141 float vector1[3]={CrossPlane[0]-CrossPlane[3],
00142 CrossPlane[1]-CrossPlane[4],
00143 CrossPlane[2]-CrossPlane[5]};
00144 float vector2[3]={CrossPlane[0]-CrossPlane[6],
00145 CrossPlane[1]-CrossPlane[7],
00146 CrossPlane[2]-CrossPlane[8]};
00147 float ortho[3]={vector1[1]*vector2[2]-vector1[2]*vector2[1],
00148 vector1[2]*vector2[0]-vector1[0]*vector2[2],
00149 vector1[0]*vector2[1]-vector1[1]*vector2[0]};
00150 float linevector[12*3]={0,0,ext[2],
00151 0,0,ext[2],
00152 0,0,ext[2],
00153 0,0,ext[2],
00154 0,ext[1],0,
00155 0,ext[1],0,
00156 0,ext[1],0,
00157 0,ext[1],0,
00158 ext[0],0,0,
00159 ext[0],0,0,
00160 ext[0],0,0,
00161 ext[0],0,0};
00162 float lineori[12*3]={ori[0],ori[1],ori[2],
00163 ori[0]+ext[0],ori[1],ori[2],
00164 ori[0]+ext[0],ori[1]+ext[1],ori[2],
00165 ori[0],ori[1]+ext[1],ori[2],
00166
00167 ori[0],ori[1],ori[2],
00168 ori[0]+ext[0],ori[1],ori[2],
00169 ori[0]+ext[0],ori[1],ori[2]+ext[2],
00170 ori[0],ori[1],ori[2]+ext[2],
00171
00172 ori[0],ori[1],ori[2],
00173 ori[0],ori[1],ori[2]+ext[2],
00174 ori[0],ori[1]+ext[1],ori[2]+ext[2],
00175 ori[0],ori[1]+ext[1],ori[2]};
00176 float polygon[4*3];
00177 int index=0;
00178
00179 for(int i=0;i<12;i++){
00180
00181 float A[9]={vector1[0],vector2[0],-linevector[3*i],
00182 vector1[1],vector2[1],-linevector[3*i+1],
00183 vector1[2],vector2[2],-linevector[3*i+2]};
00184 float b[3]={lineori[3*i]-CrossPlane[0],
00185 lineori[3*i+1]-CrossPlane[1],
00186 lineori[3*i+2]-CrossPlane[2]};
00187 float parallel=linevector[3*i]*ortho[0]
00188 +linevector[3*i+1]*ortho[1]
00189 +linevector[3*i+2]*ortho[2];
00190
00191 if(std::fabs(parallel)>tol){
00192 gauss(A,b);
00193 int k;
00194 for(int j=0;j<3;j++)
00195 polygon[index*3+j]=b[2]*linevector[3*i+j]+lineori[3*i+j];
00196
00197 if(tol<b[2] && b[2]<1.0-tol)
00198 index++;
00199 if((std::fabs(b[2])<tol || std::fabs(1.0-b[2])<tol)){
00200 if(index>0){
00201
00202 for(k=0;k<index;k++){
00203 if(std::fabs(polygon[index*3]-polygon[k*3])<tol &&
00204 std::fabs(polygon[index*3+1]-polygon[k*3+1])<tol &&
00205 std::fabs(polygon[index*3+2]-polygon[k*3+2])<tol)
00206 break;
00207 }
00208 if(k==index)
00209 index++;
00210 }
00211 else
00212 index++;
00213
00214 }
00215
00216
00217 }
00218 }
00219
00220 if(index>2){
00221 float radii = 0.;
00222 int intzero = 0;
00223 int itype = 1;
00224 int icolor = 1;
00225 int np=3*3;
00226 if(index==3)
00227 np=3;
00228 float color;
00229 switch (base::_Level){
00230 case 0:
00231 color=.6;
00232 break;
00233 case 1:
00234 color=.5;
00235 break;
00236 case 2:
00237 color=.4;
00238 break;
00239 case 3:
00240 color=.3;
00241 break;
00242 case 4:
00243 color=.2;
00244 break;
00245 default:
00246 color=.1;
00247 break;
00248 }
00249
00250 long iv[1]={1};
00251 float rv[3]={color,color,color};
00252 float saveentry[3];
00253 int opt=18;
00254 char blank=' ';
00255 V3_GETSTATE(&opt,iv,saveentry,&blank,1);
00256 V3_SETSTATE(&opt,iv,rv,&blank,1);
00257 color=0;
00258 float triang[3*3*3]={polygon[0],polygon[1],polygon[2],
00259 polygon[3],polygon[4],polygon[5],
00260 polygon[6],polygon[7],polygon[8],
00261
00262 polygon[3],polygon[4],polygon[5],
00263 polygon[6],polygon[7],polygon[8],
00264 polygon[9],polygon[10],polygon[11],
00265
00266 polygon[6],polygon[7],polygon[8],
00267 polygon[9],polygon[10],polygon[11],
00268 polygon[0],polygon[1],polygon[2]};
00269 float triangcol[3*3]={color,color,color,
00270 color,color,color,
00271 color,color,color};
00272 V3_OBJECT3D(&itype, &icolor, triang, &radii, triangcol, &intzero, &np);
00273 V3_SETSTATE(&opt,iv,saveentry,&blank,1);
00274 }
00275 }
00276
00277
00278 int FindMatchingNode(int ijk[]){
00279
00280 if (ijk[0]<base::DB().lower(0) || ijk[0]>base::DB().upper(0)+base::DB().stepsize(0) ||
00281 ijk[1]<base::DB().lower(1) || ijk[1]>base::DB().upper(1)+base::DB().stepsize(1) ||
00282 ijk[2]<base::DB().lower(2) || ijk[2]>base::DB().upper(2)+base::DB().stepsize(2))
00283 return -1;
00284
00285 if (ijk[0]%base::DB().stepsize(0) != 0 ||
00286 ijk[1]%base::DB().stepsize(1) != 0 ||
00287 ijk[2]%base::DB().stepsize(2) != 0)
00288 return -1;
00289
00290 int exti=base::DB().extents(0)+1;
00291 int extj=base::DB().extents(1)+1;
00292 return ( (ijk[2]-base::DB().lower(2))/base::DB().stepsize(2)*exti*extj +
00293 (ijk[1]-base::DB().lower(1))/base::DB().stepsize(1)*exti +
00294 (ijk[0]-base::DB().lower(0))/base::DB().stepsize(0) + 1 + base::_nodeOffset );
00295 }
00296
00297 inline const int& NodeOffset() const { return base::_nodeOffset; }
00298 inline int NodeOffset() { return base::_nodeOffset; }
00299
00300 int ExternalSurface(Coords& slc, Coords& suc, BBoxList& surflist) {
00301 for (BBox *bb = surflist.first(); bb ; bb = surflist.next())
00302 if (bb->inside(slc) && bb->inside(suc))
00303 return 0;
00304 return 1;
00305 }
00306
00307 void SetSurfaces(int dir, BBoxList& surflist, int& ksurf, int scel[][4],
00308 int set) {
00309
00310 int exti=base::DB().extents(0)+1;
00311 int extj=base::DB().extents(1)+1;
00312 int extk=base::DB().extents(2)+1;
00313 int i, j, k, kn;
00314 const Coords& base = base::DB().lower();
00315 const Coords& step = base::DB().stepsize();
00316
00317 switch (dir) {
00318 case 0:
00319 for (k = 0; k < extk-1; k++){
00320 for (j = 0; j < extj-1; j++){
00321 kn = k*exti*extj + j*exti + 1+base::_nodeOffset;
00322
00323 Coords slc(3,base(0),base(1)+ j *step(1),base(2)+ k *step(2));
00324 Coords suc(3,base(0),base(1)+(j+1)*step(1),base(2)+(k+1)*step(2));
00325 if (ExternalSurface(slc,suc,surflist)) {
00326 if (set) {
00327 scel[ksurf][0] = kn;
00328 scel[ksurf][1] = kn + exti;
00329 scel[ksurf][2] = kn + exti + exti*extj;
00330 scel[ksurf][3] = kn + exti*extj;
00331 }
00332 ksurf++;
00333 }
00334 }
00335 }
00336 break;
00337 case 1:
00338 for (k = 0; k < extk-1; k++){
00339 for (j = 0; j < extj-1; j++){
00340 kn = k*exti*extj + j*exti + exti+base::_nodeOffset;
00341
00342 Coords slc(3,base(0)+(exti-1)*step(0),base(1)+ j *step(1),
00343 base(2)+ k *step(2));
00344 Coords suc(3,base(0)+(exti-1)*step(0),base(1)+(j+1)*step(1),
00345 base(2)+(k+1)*step(2));
00346 if (ExternalSurface(slc,suc,surflist)) {
00347 if (set) {
00348 scel[ksurf][0] = kn;
00349 scel[ksurf][1] = kn + exti;
00350 scel[ksurf][2] = kn + exti + exti*extj;
00351 scel[ksurf][3] = kn + exti*extj;
00352 }
00353 ksurf++;
00354 }
00355 }
00356 }
00357 break;
00358 case 2:
00359 for (k = 0; k < extk-1; k++){
00360 for (i = 0; i < exti-1; i++){
00361 kn = k*exti*extj + i + 1+base::_nodeOffset;
00362
00363 Coords slc(3,base(0)+ i *step(0),base(1),base(2)+ k *step(2));
00364 Coords suc(3,base(0)+(i+1)*step(0),base(1),base(2)+(k+1)*step(2));
00365 if (ExternalSurface(slc,suc,surflist)) {
00366 if (set) {
00367 scel[ksurf][0] = kn;
00368 scel[ksurf][1] = kn + 1;
00369 scel[ksurf][2] = kn + 1 + exti*extj;
00370 scel[ksurf][3] = kn + exti*extj;
00371 }
00372 ksurf++;
00373 }
00374 }
00375 }
00376 break;
00377 case 3:
00378 for (k = 0; k < extk-1; k++){
00379 for (i = 0; i < exti-1; i++){
00380 kn = k*exti*extj + (extj-1)*exti + i + 1+base::_nodeOffset;
00381
00382 Coords slc(3,base(0)+ i *step(0),base(1)+(extj-1)*step(1),
00383 base(2)+ k *step(2));
00384 Coords suc(3,base(0)+(i+1)*step(0),base(1)+(extj-1)*step(1),
00385 base(2)+(k+1)*step(2));
00386 if (ExternalSurface(slc,suc,surflist)) {
00387 if (set) {
00388 scel[ksurf][0] = kn;
00389 scel[ksurf][1] = kn + 1;
00390 scel[ksurf][2] = kn + 1 + exti*extj;
00391 scel[ksurf][3] = kn + exti*extj;
00392 }
00393 ksurf++;
00394 }
00395 }
00396 }
00397 break;
00398 case 4:
00399 for (j = 0; j < extj-1; j++){
00400 for (i = 0; i < exti-1; i++){
00401 kn = j*exti + i + 1+base::_nodeOffset;
00402
00403 Coords slc(3,base(0)+ i *step(0),base(1)+ j *step(1),base(2));
00404 Coords suc(3,base(0)+(i+1)*step(0),base(1)+(j+1)*step(1),base(2));
00405 if (ExternalSurface(slc,suc,surflist)) {
00406 if (set) {
00407 scel[ksurf][0] = kn;
00408 scel[ksurf][1] = kn + 1;
00409 scel[ksurf][2] = kn + 1 + exti;
00410 scel[ksurf][3] = kn + exti;
00411 }
00412 ksurf++;
00413 }
00414 }
00415 }
00416 break;
00417 case 5:
00418 for (j = 0; j < extj-1; j++){
00419 for (i = 0; i < exti-1; i++){
00420 kn = (extk-1)*extj*exti + j*exti + i + 1+base::_nodeOffset;
00421
00422 Coords slc(3,base(0)+ i *step(0),base(1)+ j *step(1),
00423 base(2)+(extk-1)*step(2));
00424 Coords suc(3,base(0)+(i+1)*step(0),base(1)+(j+1)*step(1),
00425 base(2)+(extk-1)*step(2));
00426 if (ExternalSurface(slc,suc,surflist)) {
00427 if (set) {
00428 scel[ksurf][0] = kn;
00429 scel[ksurf][1] = kn + 1;
00430 scel[ksurf][2] = kn + 1 + exti;
00431 scel[ksurf][3] = kn + exti;
00432 }
00433 ksurf++;
00434 }
00435 }
00436 }
00437 break;
00438
00439 default:
00440 break;
00441 }
00442 }
00443
00444 void SetCrossPlane(float* cp) { CrossPlane = cp; }
00445
00446 private:
00447 void gauss(float *a,float *b) {
00448 int i,j,k;
00449 float dummy;
00450
00451 for(k=0;k<2;k++){
00452 for(i=k;i<2;i++){
00453 if(std::fabs(a[i*3])>tol)
00454 break;
00455 }
00456 for(j=0;j<3;j++){
00457 dummy=a[i*3+j];
00458 a[i*3+j]=a[k*3+j];
00459 a[k*3+j]=dummy;
00460 }
00461 dummy=b[i];
00462 b[i]=b[k];
00463 b[k]=dummy;
00464
00465 }
00466
00467 if(std::fabs(a[3])>tol){
00468 a[4]=a[3]*a[1]-a[0]*a[4];
00469 a[5]=a[3]*a[2]-a[0]*a[5];
00470 b[1]=a[3]*b[0]-a[0]*b[1];
00471 }
00472 if(std::fabs(a[6])>tol){
00473 a[7]=a[6]*a[1]-a[0]*a[7];
00474 a[8]=a[6]*a[2]-a[0]*a[8];
00475 b[2]=a[6]*b[0]-a[0]*b[2];
00476 }
00477
00478 b[2]= (a[4]*b[2]-a[7]*b[1])/(a[4]*a[8]-a[7]*a[5]);
00479 }
00480
00481 private:
00482 float* CrossPlane;
00483 };
00484
00485 #endif
