amroc/converters/Evaluator.h
Go to the documentation of this file.00001 // -*- C++ -*- 00002 00003 // Copyright (C) 2002 Ralf Deiterding 00004 // Brandenburgische Universitaet Cottbus 00005 // 00006 // Copyright (C) 2003-2007 California Institute of Technology 00007 // Ralf Deiterding, ralf@cacr.caltech.edu 00008 00009 #ifndef AMROC_EVALUATOR_H 00010 #define AMROC_EVALUATOR_H 00011 00012 #include <iostream> 00013 #include <cstdio> 00014 #include <cmath> 00015 #include <cfloat> 00016 00024 #include "EvaluatorBase.h" 00025 00033 template <class VectorType> 00034 class Evaluator : public EvaluatorBase { 00035 typedef typename VectorType::InternalDataType DataType; 00036 typedef EvaluatorBase base; 00037 00038 public: 00039 typedef GridData<VectorType,3> vector_block_type; 00040 typedef GridData<DataType,3> data_block_type; 00041 00042 Evaluator() : EvaluatorBase() { 00043 std::sprintf(title,"3D-Euler equations"); 00044 } 00045 00046 virtual void SetUp(VisualGridBase *gr) { 00047 EvaluatorBase::SetUp(gr); 00048 00049 std::sprintf(&(tkeys[ 0*LEN_TKEYS]), "Distribution "); 00050 std::sprintf(&(tkeys[ 1*LEN_TKEYS]), "Levels "); 00051 std::sprintf(&(tkeys[ 2*LEN_TKEYS]), "Density "); 00052 std::sprintf(&(tkeys[ 3*LEN_TKEYS]), "Velocity u "); 00053 std::sprintf(&(tkeys[ 4*LEN_TKEYS]), "Velocity v "); 00054 std::sprintf(&(tkeys[ 5*LEN_TKEYS]), "Velocity w "); 00055 std::sprintf(&(tkeys[ 6*LEN_TKEYS]), "Total Energy "); 00056 std::sprintf(&(tkeys[ 7*LEN_TKEYS]), "|Velocity| "); 00057 std::sprintf(&(tkeys[ 8*LEN_TKEYS]), "Flow Vectors "); 00058 std::sprintf(&(tkeys[ 9*LEN_TKEYS]), "Schl.-Plot Density "); 00059 std::sprintf(&(tkeys[10*LEN_TKEYS]), "Schl.-Plot Velocity u "); 00060 std::sprintf(&(tkeys[11*LEN_TKEYS]), "Schl.-Plot Velocity v "); 00061 std::sprintf(&(tkeys[12*LEN_TKEYS]), "Schl.-Plot Velocity w "); 00062 std::sprintf(&(tkeys[13*LEN_TKEYS]), "Schl.-Plot Total Energy "); 00063 std::sprintf(&(tkeys[14*LEN_TKEYS]), "Schl.-Plot |Velocity| "); 00064 std::sprintf(&(tkeys[15*LEN_TKEYS]), "Vorticity of w and v "); 00065 std::sprintf(&(tkeys[16*LEN_TKEYS]), "Vorticity of u and w "); 00066 std::sprintf(&(tkeys[17*LEN_TKEYS]), "Vorticity of v and u "); 00067 00068 base::fkeys[0] = base::Grid->FKeys(); 00069 base::fkeys[1] = base::Grid->FKeys(); 00070 base::fkeys[2] = base::Grid->FKeys(); 00071 base::fkeys[3] = base::Grid->FKeys(); 00072 base::fkeys[4] = base::Grid->FKeys(); 00073 base::fkeys[5] = base::Grid->FKeys(); 00074 base::fkeys[6] = base::Grid->FKeys(); 00075 base::fkeys[7] = base::Grid->FKeys(); 00076 base::fkeys[8] = 2; 00077 base::fkeys[9] = base::Grid->FKeys(); 00078 base::fkeys[10] = base::Grid->FKeys(); 00079 base::fkeys[11] = base::Grid->FKeys(); 00080 base::fkeys[12] = base::Grid->FKeys(); 00081 base::fkeys[13] = base::Grid->FKeys(); 00082 base::fkeys[14] = base::Grid->FKeys(); 00083 base::fkeys[15] = base::Grid->FKeys(); 00084 base::fkeys[16] = base::Grid->FKeys(); 00085 base::fkeys[17] = base::Grid->FKeys(); 00086 00087 base::ikeys[0] = 105; base::ikeys[1] = 115; base::ikeys[2] = 100; base::ikeys[3] = 117; 00088 base::ikeys[4] = 118; base::ikeys[5] = 119; base::ikeys[6] = 101; base::ikeys[7] = 108; 00089 base::ikeys[8] = 102; base::ikeys[9] = 68; base::ikeys[10] = 85; base::ikeys[11] = 86; 00090 base::ikeys[12] = 87; base::ikeys[13] = 69; base::ikeys[14] = 76; base::ikeys[15] = 120; 00091 base::ikeys[16] = 121; base::ikeys[17] = 122; 00092 } 00093 00094 virtual int NKeys() const { return 18; } 00095 00096 virtual void SetScalCells(int key,float v[],int& offset, 00097 vector_block_type& DB, float* dx, bool* CompRead) { 00098 int idx=0; 00099 int keym1 = key-1; 00100 int ac; 00101 BeginFastIndex3(dat, DB.bbox(), DB.data(), VectorType); 00102 switch (key) { 00103 case 3: 00104 case 4: 00105 case 5: 00106 case 6: 00107 case 7: 00108 case 10: 00109 case 11: 00110 case 12: 00111 case 13: 00112 case 14: 00113 ac = ((key>=10 && key<=14) ? 7 : 0); 00114 if (!CompRead[key-(3+ac)]) 00115 return; 00116 for_3 (i,j,k,DB.bbox(),DB.bbox().stepsize()) 00117 v[idx+offset]= FastIndex3(dat,i,j,k)(key-(3+ac)); 00118 mvals[2*keym1] = Min(mvals[2*keym1], v[idx+offset]); 00119 mvals[2*keym1+1] = Max(mvals[2*keym1+1], v[idx+offset]); 00120 idx++; 00121 end_for 00122 break; 00123 00124 case 8: 00125 case 15: 00126 case 16: 00127 case 17: 00128 case 18: 00129 float u1, u2, u3; 00130 for_3 (i,j,k,DB.bbox(),DB.bbox().stepsize()) 00131 u1=FastIndex3(dat,i,j,k)(1); 00132 u2=FastIndex3(dat,i,j,k)(2); 00133 u3=FastIndex3(dat,i,j,k)(3); 00134 v[idx+offset]=std::sqrt(u1*u1+u2*u2+u3*u3); 00135 mvals[2*keym1] = Min(mvals[2*keym1], v[idx+offset]); 00136 mvals[2*keym1+1] = Max(mvals[2*keym1+1], v[idx+offset]); 00137 idx++; 00138 end_for 00139 break; 00140 00141 default: 00142 break; 00143 } 00144 EndFastIndex3(dat); 00145 offset+=idx; 00146 00147 if ((key>=10 && key<=14) || (key>=15 && key<=18)) 00148 std::cerr << &(tkeys[keym1*LEN_TKEYS]) 00149 << " is NOT supported for Type=6. Use Type=1 instead!" 00150 << std::endl; 00151 } 00152 00153 virtual void SetScalNodes(int key,float v[],int& offset, 00154 vector_block_type& DB, const BBox& bbox, 00155 float* dx, bool* CompRead) { 00156 int idx=0; 00157 int keym1 = key-1; 00158 BeginFastIndex3(dat, DB.bbox(), DB.data(), VectorType); 00159 const Coords& step = bbox.stepsize(); 00160 switch (key) { 00161 case 3: 00162 case 4: 00163 case 5: 00164 case 6: 00165 case 7: 00166 { 00167 if (!CompRead[key-3]) 00168 return; 00169 for_3 (i,j,k,bbox,step) 00170 BBox AddOver(3,i-step(0),j-step(1),k-step(2), 00171 i,j,k,step(0),step(1),step(2)); 00172 AddOver *= DB.bbox(); 00173 00174 float c=0.; 00175 v[idx+offset] = 0.; 00176 for_3 (l, m, n, AddOver, step) 00177 if (FastIndex3(dat,l,m,n)(key-3) < FLT_MAX) { 00178 v[idx+offset] += FastIndex3(dat,l,m,n)(key-3); 00179 c += 1.; 00180 } 00181 end_for 00182 if (c>0) v[idx+offset] /= c; 00183 else v[idx+offset] = FLT_MAX; 00184 00185 mvals[2*keym1] = Min(mvals[2*keym1], v[idx+offset]); 00186 mvals[2*keym1+1] = Max(mvals[2*keym1+1], v[idx+offset]); 00187 idx++; 00188 end_for 00189 break; 00190 } 00191 case 10: 00192 case 11: 00193 case 12: 00194 case 13: 00195 case 14: 00196 { 00197 if (!CompRead[key-10]) 00198 return; 00199 for_3 (i,j,k,bbox,step) 00200 BBox AddOver(3,i-step(0),j-step(1),k-step(2), 00201 i,j,k,step(0),step(1),step(2)); 00202 AddOver *= DB.bbox(); 00203 00204 data_block_type DBHelp(AddOver); 00205 BeginFastIndex3(help, DBHelp.bbox(), DBHelp.data(), DataType); 00206 for_3 (l, m, n, AddOver, step) 00207 FastIndex3(help,l,m,n) = FastIndex3(dat,l,m,n)(key-10); 00208 end_for 00209 EndFastIndex3(help); 00210 v[idx+offset] = SetScalGradNodes(DBHelp,dx); 00211 00212 mvals[2*keym1] = Min(mvals[2*keym1], v[idx+offset]); 00213 mvals[2*keym1+1] = Max(mvals[2*keym1+1], v[idx+offset]); 00214 idx++; 00215 end_for 00216 break; 00217 } 00218 case 8: 00219 { 00220 float u1, u2, u3; 00221 for_3 (i, j, k, bbox, step) 00222 BBox AddOver(3,i-step(0),j-step(1),k-step(2), 00223 i,j,k,step(0),step(1),step(2)); 00224 AddOver *= DB.bbox(); 00225 00226 float c=0.; 00227 v[idx+offset] = 0.; 00228 for_3 (l, m, n, AddOver, step) 00229 u1 = FastIndex3(dat,l,m,n)(1); 00230 u2 = FastIndex3(dat,l,m,n)(2); 00231 u3 = FastIndex3(dat,l,m,n)(3); 00232 if (u1<FLT_MAX && u2<FLT_MAX && u3<FLT_MAX) { 00233 v[idx+offset] += std::sqrt(u1*u1+u2*u2+u3*u3); 00234 c += 1.; 00235 } 00236 end_for 00237 if (c>0) v[idx+offset] /= c; 00238 else v[idx+offset] = FLT_MAX; 00239 00240 mvals[2*keym1] = Min(mvals[2*keym1], v[idx+offset]); 00241 mvals[2*keym1+1] = Max(mvals[2*keym1+1], v[idx+offset]); 00242 idx++; 00243 end_for 00244 break; 00245 } 00246 case 15: 00247 { 00248 float u1, u2, u3; 00249 for_3 (i, j, k, bbox, step) 00250 BBox AddOver(3,i-step(0),j-step(1),k-step(2), 00251 i,j,k,step(0),step(1),step(2)); 00252 AddOver *= DB.bbox(); 00253 00254 data_block_type DBHelp(AddOver); 00255 DBHelp = FLT_MAX; 00256 BeginFastIndex3(help, DBHelp.bbox(), DBHelp.data(), DataType); 00257 for_3 (l, m, n, AddOver, step) 00258 u1 = FastIndex3(dat,l,m,n)(1); 00259 u2 = FastIndex3(dat,l,m,n)(2); 00260 u3 = FastIndex3(dat,l,m,n)(3); 00261 if (u1<FLT_MAX && u2<FLT_MAX && u3<FLT_MAX) 00262 FastIndex3(help,l,m,n) = std::sqrt(u1*u1+u2*u2+u3*u3); 00263 end_for 00264 EndFastIndex3(help); 00265 v[idx+offset] = SetScalGradNodes(DBHelp,dx); 00266 00267 mvals[2*keym1] = Min(mvals[2*keym1], v[idx+offset]); 00268 mvals[2*keym1+1] = Max(mvals[2*keym1+1], v[idx+offset]); 00269 idx++; 00270 end_for 00271 break; 00272 } 00273 case 16: 00274 { 00275 if (!CompRead[2] || !CompRead[3]) 00276 return; 00277 for_3 (i,j,k,bbox,step) 00278 BBox AddOver(3,i-step(0),j-step(1),k-step(2), 00279 i,j,k,step(0),step(1),step(2)); 00280 AddOver *= DB.bbox(); 00281 00282 data_block_type u3(AddOver), u2(AddOver); 00283 BeginFastIndex3(u2, u2.bbox(), u2.data(), DataType); 00284 BeginFastIndex3(u3, u3.bbox(), u3.data(), DataType); 00285 for_3 (l, m, n, AddOver, step) 00286 FastIndex3(u2,l,m,n) = FastIndex3(dat,l,m,n)(2); 00287 FastIndex3(u3,l,m,n) = FastIndex3(dat,l,m,n)(3); 00288 end_for 00289 EndFastIndex3(u2); 00290 EndFastIndex3(u3); 00291 v[idx+offset] = ScalDerivNodes(u3,dx,1)-ScalDerivNodes(u2,dx,2); 00292 00293 mvals[2*keym1] = Min(mvals[2*keym1], v[idx+offset]); 00294 mvals[2*keym1+1] = Max(mvals[2*keym1+1], v[idx+offset]); 00295 idx++; 00296 end_for 00297 break; 00298 } 00299 case 17: 00300 { 00301 if (!CompRead[1] || !CompRead[3]) 00302 return; 00303 for_3 (i,j,k,bbox,step) 00304 BBox AddOver(3,i-step(0),j-step(1),k-step(2), 00305 i,j,k,step(0),step(1),step(2)); 00306 AddOver *= DB.bbox(); 00307 00308 data_block_type u1(AddOver), u3(AddOver); 00309 BeginFastIndex3(u1, u1.bbox(), u1.data(), DataType); 00310 BeginFastIndex3(u3, u3.bbox(), u3.data(), DataType); 00311 for_3 (l, m, n, AddOver, step) 00312 FastIndex3(u1,l,m,n) = FastIndex3(dat,l,m,n)(1); 00313 FastIndex3(u3,l,m,n) = FastIndex3(dat,l,m,n)(3); 00314 end_for 00315 EndFastIndex3(u1); 00316 EndFastIndex3(u3); 00317 v[idx+offset] = ScalDerivNodes(u1,dx,2)-ScalDerivNodes(u3,dx,0); 00318 00319 mvals[2*keym1] = Min(mvals[2*keym1], v[idx+offset]); 00320 mvals[2*keym1+1] = Max(mvals[2*keym1+1], v[idx+offset]); 00321 idx++; 00322 end_for 00323 break; 00324 } 00325 case 18: 00326 { 00327 if (!CompRead[1] || !CompRead[2]) 00328 return; 00329 for_3 (i,j,k,bbox,step) 00330 BBox AddOver(3,i-step(0),j-step(1),k-step(2), 00331 i,j,k,step(0),step(1),step(2)); 00332 AddOver *= DB.bbox(); 00333 00334 data_block_type u1(AddOver), u2(AddOver); 00335 BeginFastIndex3(u1, u1.bbox(), u1.data(), DataType); 00336 BeginFastIndex3(u2, u2.bbox(), u2.data(), DataType); 00337 for_3 (l, m, n, AddOver, step) 00338 FastIndex3(u1,l,m,n) = FastIndex3(dat,l,m,n)(1); 00339 FastIndex3(u2,l,m,n) = FastIndex3(dat,l,m,n)(2); 00340 end_for 00341 EndFastIndex3(u1); 00342 EndFastIndex3(u2); 00343 v[idx+offset] = ScalDerivNodes(u2,dx,0)-ScalDerivNodes(u1,dx,1); 00344 00345 mvals[2*keym1] = Min(mvals[2*keym1], v[idx+offset]); 00346 mvals[2*keym1+1] = Max(mvals[2*keym1+1], v[idx+offset]); 00347 idx++; 00348 end_for 00349 break; 00350 } 00351 00352 default: 00353 break; 00354 } 00355 EndFastIndex3(dat); 00356 offset+=idx; 00357 } 00358 00359 protected: 00360 DataType ScalDerivNodes(data_block_type& DB, float* dx, const int d) { 00361 float grad = 0.0; 00362 float c=0.0; 00363 BeginFastIndex3(dat, DB.bbox(), DB.data(), DataType); 00364 const BBox& AddOver = DB.bbox(); 00365 float wdx[3]; 00366 float v1, v2; 00367 for (int i=0; i<3; i++) 00368 wdx[i] = dx[i]*AddOver.stepsize(i); 00369 00370 switch(d) { 00371 case 0: 00372 if (AddOver.extents(0) == 2) { 00373 for (int m=AddOver.lower(1); m<=AddOver.upper(1); 00374 m+=AddOver.stepsize(1)) { 00375 for (int n=AddOver.lower(2); n<=AddOver.upper(2); 00376 n+=AddOver.stepsize(2)) { 00377 v1 = FastIndex3(dat,AddOver.lower(0),m,n); 00378 v2 = FastIndex3(dat,AddOver.upper(0),m,n); 00379 if (v1<FLT_MAX && v2<FLT_MAX) { 00380 grad += (v2-v1)/wdx[0]; 00381 c += 1.0; 00382 } 00383 else if (v1<FLT_MAX || v2<FLT_MAX) c += 1.0; 00384 } 00385 } 00386 if (c>0.0) grad /= c; 00387 else grad = FLT_MAX; 00388 } 00389 break; 00390 case 1: 00391 if (AddOver.extents(1) == 2) { 00392 for (int l=AddOver.lower(0); l<=AddOver.upper(0); 00393 l+=AddOver.stepsize(0)) { 00394 for (int n=AddOver.lower(2); n<=AddOver.upper(2); 00395 n+=AddOver.stepsize(2)) { 00396 v1 = FastIndex3(dat,l,AddOver.lower(1),n); 00397 v2 = FastIndex3(dat,l,AddOver.upper(1),n); 00398 if (v1<FLT_MAX && v2<FLT_MAX) { 00399 grad += (v2-v1)/wdx[1]; 00400 c += 1.0; 00401 } 00402 else if (v1<FLT_MAX || v2<FLT_MAX) c += 1.0; 00403 } 00404 } 00405 if (c>0.0) grad /= c; 00406 else grad = FLT_MAX; 00407 } 00408 break; 00409 case 2: 00410 if (AddOver.extents(2) == 2) { 00411 for (int l=AddOver.lower(0); l<=AddOver.upper(0); 00412 l+=AddOver.stepsize(0)) { 00413 for (int m=AddOver.lower(1); m<=AddOver.upper(1); 00414 m+=AddOver.stepsize(1)) { 00415 v1 = FastIndex3(dat,l,m,AddOver.lower(2)); 00416 v2 = FastIndex3(dat,l,m,AddOver.upper(2)); 00417 if (v1<FLT_MAX && v2<FLT_MAX) { 00418 grad += (v2-v1)/wdx[2]; 00419 c += 1.0; 00420 } 00421 else if (v1<FLT_MAX || v2<FLT_MAX) c += 1.0; 00422 } 00423 } 00424 if (c>0.0) grad /= c; 00425 else grad = FLT_MAX; 00426 } 00427 break; 00428 default: 00429 break; 00430 } 00431 EndFastIndex3(dat); 00432 00433 return grad; 00434 } 00435 00436 DataType SetScalGradNodes(data_block_type& DB, float* dx) { 00437 float xgrad, ygrad, zgrad; 00438 xgrad = ScalDerivNodes(DB,dx,0); 00439 ygrad = ScalDerivNodes(DB,dx,1); 00440 zgrad = ScalDerivNodes(DB,dx,2); 00441 if (xgrad<FLT_MAX && ygrad<FLT_MAX && zgrad<FLT_MAX) 00442 return std::sqrt(xgrad*xgrad+ygrad*ygrad+zgrad*zgrad); 00443 else 00444 return FLT_MAX; 00445 } 00446 }; 00447 00448 00449 #endif
Generated on Thu Jun 30 02:19:02 2016 for AMROC's HDF Tools by
1.6.3