00001
00002
00003
00004
00005
00006
00007
00008
00009 #ifndef AMROC_MULTICOMP_EULER_EVALUATOR_H
00010 #define AMROC_MULTICOMP_EULER_EVALUATOR_H
00011
00019 #include <iostream>
00020 #include <string>
00021 #include <cstdio>
00022 #include <cmath>
00023 #include <cfloat>
00024
00025 #include "Evaluator.h"
00026
00033 template <class VectorType>
00034 class MultiCompEulerEvaluator : public Evaluator<VectorType> {
00035 typedef Evaluator<VectorType> base;
00036 typedef typename VectorType::InternalDataType DataType;
00037
00038 public:
00039 typedef typename base::vector_block_type vector_block_type;
00040 typedef typename base::data_block_type data_block_type;
00041
00042 MultiCompEulerEvaluator() : base(), RU(1.), PA(1.) {
00043 std::sprintf(base::title,"3D-Euler equations - Multiple Components");
00044 VectorType dummy;
00045 Nsp = dummy.length()-8;
00046 SpName = new std::string[Nsp];
00047 Wk = new DataType[Nsp];
00048 for (register int i=0; i<Nsp; i++) Wk[i] = 1.;
00049 }
00050
00051 ~MultiCompEulerEvaluator() {
00052 if (!SpName) delete [] SpName;
00053 if (!Wk) delete [] Wk;
00054 }
00055
00056 virtual void update() {
00057 ControlDevice ChemCtrl(GetFileControlDevice("chem.dat",""));
00058 RegisterAt(ChemCtrl,"RU",RU);
00059 RegisterAt(ChemCtrl,"PA",PA);
00060 for (int k=0; k<Nsp; k++) {
00061 char RegName[16];
00062 std::sprintf(RegName,"Sp(%02d)",k+1);
00063 RegisterAt(ChemCtrl,RegName,SpName[k]);
00064 std::sprintf(RegName,"W(%02d)",k+1);
00065 RegisterAt(ChemCtrl,RegName,Wk[k]);
00066 }
00067 ChemCtrl.update();
00068 }
00069
00070 virtual void SetUp(VisualGridBase *gr) {
00071 base::SetUp(gr);
00072 const int bk=base::NKeys();
00073 std::sprintf(&(base::tkeys[ bk *LEN_TKEYS]), "Temperature ");
00074 std::sprintf(&(base::tkeys[(bk+1)*LEN_TKEYS]), "Pressure ");
00075 std::sprintf(&(base::tkeys[(bk+2)*LEN_TKEYS]), "Gamma ");
00076 std::sprintf(&(base::tkeys[(bk+3)*LEN_TKEYS]), "Speed of Sound ");
00077 std::sprintf(&(base::tkeys[(bk+4)*LEN_TKEYS]), "Schl.-Plot Temperature ");
00078 std::sprintf(&(base::tkeys[(bk+5)*LEN_TKEYS]), "Schl.-Plot Pressure ");
00079 std::sprintf(&(base::tkeys[(bk+6)*LEN_TKEYS]), "Schl.-Plot Gamma ");
00080 std::sprintf(&(base::tkeys[(bk+7)*LEN_TKEYS]), "Schl.-Plot Sound-Speed ");
00081 base::fkeys[bk ] = base::Grid->FKeys();
00082 base::fkeys[bk+1] = base::Grid->FKeys();
00083 base::fkeys[bk+2] = base::Grid->FKeys();
00084 base::fkeys[bk+3] = base::Grid->FKeys();
00085 base::fkeys[bk+4] = base::Grid->FKeys();
00086 base::fkeys[bk+5] = base::Grid->FKeys();
00087 base::fkeys[bk+6] = base::Grid->FKeys();
00088 base::fkeys[bk+7] = base::Grid->FKeys();
00089 base::ikeys[bk ] = 116;
00090 base::ikeys[bk+1] = 112;
00091 base::ikeys[bk+2] = 103;
00092 base::ikeys[bk+3] = 99;
00093 base::ikeys[bk+4] = 84;
00094 base::ikeys[bk+5] = 80;
00095 base::ikeys[bk+6] = 71;
00096 base::ikeys[bk+7] = 67;
00097
00098 for (int i=0; i<Nsp; i++) {
00099 std::sprintf(&(base::tkeys[(bk+8+i)*LEN_TKEYS]),"Y_%02d - %-17s",i+1,
00100 SpName[i].c_str());
00101 base::fkeys[bk+8+i] = base::Grid->FKeys();
00102 base::ikeys[bk+8+i] = 48+i;
00103 }
00104 }
00105
00106 virtual int NKeys() const { return base::NKeys()+8+Nsp; }
00107
00108 virtual void SetScalCells(int key,float v[],int& offset,
00109 vector_block_type& DB, float* dx, bool* CompRead) {
00110 int idx=0;
00111 int keym1 = key-1;
00112 const int bk=base::NKeys();
00113 BeginFastIndex3(dat, DB.bbox(), DB.data(), VectorType);
00114
00115 int ac=-1;
00116 if (key>=bk+1 && key<=bk+3)
00117 ac=key-base::NKeys()-1+5;
00118 if (key>=bk+5 && key<=bk+7)
00119 ac=key-base::NKeys()-1+5-4;
00120 else if(key>=bk+8 && key<=NKeys())
00121 ac=key-base::NKeys()-9+8;
00122
00123 if (ac>=0) {
00124 if (!CompRead[ac])
00125 return;
00126 for_3 (i,j,k,DB.bbox(),DB.bbox().stepsize())
00127 v[idx+offset]= FastIndex3(dat,i,j,k)(ac);
00128 base::mvals[2*keym1] = Min(base::mvals[2*keym1], v[idx+offset]);
00129 base::mvals[2*keym1+1] = Max(base::mvals[2*keym1+1], v[idx+offset]);
00130 idx++;
00131 end_for
00132 }
00133 else if (key==bk+4 || key==bk+8) {
00134 if (!CompRead[0] || !CompRead[6] || !CompRead[7])
00135 return;
00136
00137 DataType rho,p,gamma;
00138 for_3 (i,j,k,DB.bbox(),DB.bbox().stepsize())
00139 rho =FastIndex3(dat,i,j,k)(0);
00140 p =FastIndex3(dat,i,j,k)(6);
00141 gamma=FastIndex3(dat,i,j,k)(7);
00142 v[idx+offset]=std::sqrt(gamma*p/rho);
00143 base::mvals[2*keym1] = Min(base::mvals[2*keym1], v[idx+offset]);
00144 base::mvals[2*keym1+1] = Max(base::mvals[2*keym1+1], v[idx+offset]);
00145 idx++;
00146 end_for
00147 }
00148 else if (key<=bk)
00149 base::SetScalCells(key,v,offset,DB,dx,CompRead);
00150
00151 EndFastIndex3(dat);
00152 offset+=idx;
00153
00154 if (key>=bk+5 && key<=bk+8)
00155 std::cerr << &(base::tkeys[keym1*LEN_TKEYS])
00156 << " is NOT supported for Type=6. Use Type=1 instead!"
00157 << std::endl;
00158 }
00159
00160
00161
00162 virtual void SetScalNodes(int key,float v[],int& offset,
00163 vector_block_type& DB, const BBox& bbox,
00164 float* dx, bool* CompRead) {
00165 int idx=0;
00166 int keym1 = key-1;
00167 const int bk=base::NKeys();
00168 BeginFastIndex3(dat, DB.bbox(), DB.data(), VectorType);
00169 const Coords& step = bbox.stepsize();
00170
00171 int ac=-1;
00172 if (key>=bk+1 && key<=bk+3)
00173 ac=key-base::NKeys()-1+5;
00174 else if(key>=bk+8 && key<=NKeys())
00175 ac=key-base::NKeys()-9+8;
00176
00177 if (ac>=0) {
00178 if (!CompRead[ac])
00179 return;
00180 for_3 (i,j,k,bbox,step)
00181 BBox AddOver(3,i-step(0),j-step(1),k-step(2),
00182 i,j,k,step(0),step(1),step(2));
00183 AddOver *= DB.bbox();
00184
00185 float c=0;
00186 v[idx+offset] = 0.0;
00187 for_3 (l, m, n, AddOver, step)
00188 if (FastIndex3(dat,l,m,n)(ac)<FLT_MAX) {
00189 v[idx+offset] += FastIndex3(dat,l,m,n)(ac);
00190 c += 1.0;
00191 }
00192 end_for
00193 if (c>0) v[idx+offset] /= c;
00194
00195 base::mvals[2*keym1] = Min(base::mvals[2*keym1], v[idx+offset]);
00196 base::mvals[2*keym1+1] = Max(base::mvals[2*keym1+1], v[idx+offset]);
00197 idx++;
00198 end_for
00199 }
00200
00201 if (key==bk+4) {
00202 if (!CompRead[0] || !CompRead[6] || !CompRead[7])
00203 return;
00204
00205 DataType rho,p,gamma;
00206 for_3 (i,j,k,bbox,step)
00207 BBox AddOver(3,i-step(0),j-step(1),k-step(2),
00208 i,j,k,step(0),step(1),step(2));
00209 AddOver *= DB.bbox();
00210
00211 float c=0;
00212 v[idx+offset] = 0.0;
00213 for_3 (l, m, n, AddOver, step)
00214 rho =FastIndex3(dat,l,m,n)(0);
00215 p =FastIndex3(dat,l,m,n)(6);
00216 gamma=FastIndex3(dat,l,m,n)(7);
00217 if (rho<FLT_MAX && p<FLT_MAX && gamma<FLT_MAX) {
00218 v[idx+offset] += std::sqrt(gamma*p/rho);
00219 c += 1.0;
00220 }
00221 end_for
00222 if (c>0) v[idx+offset] /= c;
00223
00224 base::mvals[2*keym1] = Min(base::mvals[2*keym1], v[idx+offset]);
00225 base::mvals[2*keym1+1] = Max(base::mvals[2*keym1+1], v[idx+offset]);
00226 idx++;
00227 end_for
00228 }
00229
00230 if (key>=bk+5 && key<=bk+7) {
00231 ac=key-base::NKeys()-5+5;
00232
00233 if (!CompRead[ac])
00234 return;
00235 for_3 (i,j,k,bbox,step)
00236 BBox AddOver(3,i-step(0),j-step(1),k-step(2),
00237 i,j,k,step(0),step(1),step(2));
00238 AddOver *= DB.bbox();
00239
00240 data_block_type DBHelp(AddOver);
00241 BeginFastIndex3(help, DBHelp.bbox(), DBHelp.data(), DataType);
00242 for_3 (l, m, n, AddOver, step)
00243 FastIndex3(help,l,m,n) = FastIndex3(dat,l,m,n)(ac);
00244 end_for
00245 EndFastIndex3(help);
00246 v[idx+offset] = base::SetScalGradNodes(DBHelp,dx);
00247
00248 base::mvals[2*keym1] = Min(base::mvals[2*keym1], v[idx+offset]);
00249 base::mvals[2*keym1+1] = Max(base::mvals[2*keym1+1], v[idx+offset]);
00250 idx++;
00251 end_for
00252 }
00253
00254 if (key==bk+8) {
00255 if (!CompRead[0] || !CompRead[6] || !CompRead[7])
00256 return;
00257
00258 DataType rho,p,gamma;
00259 for_3 (i,j,k,bbox,step)
00260 BBox AddOver(3,i-step(0),j-step(1),k-step(2),
00261 i,j,k,step(0),step(1),step(2));
00262 AddOver *= DB.bbox();
00263
00264 data_block_type DBHelp(AddOver);
00265 DBHelp = FLT_MAX;
00266 BeginFastIndex3(help, DBHelp.bbox(), DBHelp.data(), DataType);
00267 for_3 (l, m, n, AddOver, step)
00268 rho =FastIndex3(dat,l,m,n)(0);
00269 p =FastIndex3(dat,l,m,n)(6);
00270 gamma=FastIndex3(dat,l,m,n)(7);
00271 if (rho<FLT_MAX && p<FLT_MAX && gamma<FLT_MAX)
00272 FastIndex3(help,l,m,n) = std::sqrt(gamma*p/rho);
00273 end_for
00274 EndFastIndex3(help);
00275 v[idx+offset] = base::SetScalGradNodes(DBHelp,dx);
00276
00277 base::mvals[2*keym1] = Min(base::mvals[2*keym1], v[idx+offset]);
00278 base::mvals[2*keym1+1] = Max(base::mvals[2*keym1+1], v[idx+offset]);
00279 idx++;
00280 end_for
00281 }
00282
00283 if (key<=bk)
00284 base::SetScalNodes(key,v,offset,DB,bbox,dx,CompRead);
00285
00286 EndFastIndex3(dat);
00287 offset+=idx;
00288 }
00289
00290 private:
00291 int Nsp;
00292 DataType RU, PA;
00293 std::string* SpName;
00294 DataType* Wk;
00295 };
00296
00297 #endif
