SparseToOmegaVTK.cc

00001 #include<stdlib.h>
00002 #include<stdio.h>
00003 #include"Math.hpp"
00004 #include<string.h>
00005 #include<sys/types.h>
00006 #include<sys/stat.h>
00007 
00008 #include "AdaptiveData.hpp"
00009 #include "UniformData.hpp"
00010 
00011 int debugRefine ;
00012 
00013 # define D 2
00014 
00015 int main(int argc,char **argv) {
00016   int i,j,k,agc=1 ;
00017 
00018   char ufile[D][256],ofile[256],cmd[1000] ;
00019   int  L[3]={0},A[3],E[3],a[3]={0},e[3]={0},N=6,smooth_type=-1,component=-1,statdir=-1 ;
00020   bool smooth=false, restrict=false, onefile=false, absomega=true, writeomega=false, stat=false ;
00021 
00022   // parse args
00023   if (argc<6) {
00024     printf("SparseToOmegaVTK (%dD)\n",D) ;
00025     printf("      -L <L0> <L1> %s   // refinement level for full grid on which the rotation is calculated\n",(D==2) ? "" : "<L2>") ;
00026     printf("      {-r <a0> <a1> %s <e0> <e1> %s }           // (optional) just write the values of a subset of the above grid\n", (D==2) ? "" : "<a2>" , (D==2) ? "" : "<e2>") ;
00027     printf("      {-s <smooth_typ>}   default:0 (no smooth) // (optional) additional smoothing of nodal values\n") ;
00028     printf("      {-om <component>}            abs|Omega|   // (optional) compute/write which component (0..%d) of the rotation\n",D-1) ;
00029     printf("                                                //  if -om is not specified, |abs(Omega)| is computed\n") ;
00030     printf("      {-N <Interpolet-type>}  default: 6        // which wavelets have been used for the simulations\n") ;
00031     printf("      -u3 <fileu0> <fileu1> %s  |   -u  <fileu> // read velocity from this(these) file(s)\n", (D==2) ? "" : "<file2>") ;
00032     printf("      {-o <outfile>}                            // write rotation or its components to this file\n") ;
00033     printf("      {-stat <dir>} just write the correlation mean(Omega_i Omega_j). The averages are computed over slices perpendicular to <dir>\n") ;
00034     exit(-1) ;
00035   }
00036 
00037   while (agc<argc) {
00038     sscanf(argv[agc++],"%s",cmd) ;
00039     if (strcmp(cmd,"-L")==0) { // read L , calculation grid
00040       for (i=0; i<D; i++) sscanf(argv[agc++],"%d",&L[i]) ;
00041       if (restrict==false) for (i=0; i<D; i++) { a[i]=0; e[i]=1<<L[i] ;}  
00042     }
00043 
00044     if (strcmp(cmd,"-r")==0) { // read restriction       
00045       for (i=0; i<D; i++) {
00046         sscanf(argv[agc++],"%d",&a[i]) ;
00047         a[i]=max(a[i],0) ;
00048       }
00049       for (i=0; i<D; i++) {
00050         sscanf(argv[agc++],"%d",&e[i]) ;
00051         e[i]=min(e[i],1<<L[i]) ;
00052       }
00053       restrict=true ;
00054     }
00055 
00056     if (strcmp(cmd,"-s")==0) { // read smooth type
00057       sscanf(argv[agc++],"%d",&smooth_type) ;
00058       smooth=true ;
00059       smooth_type=max(0,min(smooth_type,0)) ;
00060     }
00061 
00062     if (strcmp(cmd,"-om")==0) { // read component
00063       sscanf(argv[agc++],"%d",&component) ;
00064       absomega=false;
00065       component=max(0,min(component,2)) ;
00066     }
00067 
00068     if (strcmp(cmd,"-N")==0) { // read Interpolet type
00069       sscanf(argv[agc++],"%d",&N) ;
00070       N= N & 0xfe ; //make even
00071       N=max(2,min(N,6)) ;
00072     }
00073    
00074     if (strcmp(cmd,"-u")==0) { // read velocity file
00075       sscanf(argv[agc++],"%s",ufile[0]) ;
00076       onefile=true ;
00077     }    
00078 
00079     if (strcmp(cmd,"-u3")==0) { // read velocity file
00080       for (i=0; i<D; i++) sscanf(argv[agc++],"%s",ufile[i]) ;
00081        onefile=false ;
00082     }
00083 
00084     if (strcmp(cmd,"-o")==0) { // read out file
00085       sscanf(argv[agc++],"%s",ofile) ;
00086       writeomega=true ;
00087     }
00088 
00089     if (strcmp(cmd,"-stat")==0) { // read out file
00090       sscanf(argv[agc++],"%d",&statdir) ;
00091       stat=true ;
00092     }
00093 
00094   }
00095 
00096   // print out parsed args
00097   printf("read velocity from file(s): ") ;
00098   for (i=0; i< (onefile?1 : D); i++) printf("%s ",ufile[i]) ;
00099   printf("\n") ;
00100   //
00101   printf("write result in %s\n",ofile) ;
00102   printf("for calculations/IO use finest mesh with levels L[3]= %d %d %d\n",L[0],L[1],L[2]) ;
00103   //
00104   if (restrict)
00105     printf("just write data of [%d,%d]x[%d,%d]x[%d,%d]\n",a[0],e[0],a[1],e[1],a[2],e[2]) ; 
00106   else 
00107     printf("write data of complete mesh\n") ;
00108   //
00109   if (smooth)
00110     printf("smooth data using type %d\n",smooth_type) ;
00111   else 
00112     printf("no smooth\n") ;
00113   //
00114   if (absomega)
00115     printf("write |Omega|\n") ;
00116   else 
00117     printf("write Omega_%d\n",component) ;
00118  
00119   if (stat) printf ("just write statistics perpendicular to %d\n",statdir) ;
00120   
00121   printf("use Interpolets %d\n",N) ;
00122 
00123 
00124   // GO
00125   Wavelets WC("Interpolet",N) ; 
00126   int BC[D][2] ;
00127   for (i=0;i<D; i++) { BC[i][0]=-1 ; BC[i][1]=PERIODIC ;} 
00128 
00129   for (i=0; i<D; i++) {A[i]=0; E[i]=1<<L[i] ;}
00130   UniformData<D> M(A,E,&WC) , W(A,E,&WC) ;
00131   M.SetBoundaryConditions(BC) ;
00132   W.SetBoundaryConditions(BC) ;
00133   AdaptiveGrid<D> G(&WC,false) ;
00134   G.SetPeriodicConditions(BC) ;
00135 
00136   AdaptiveData<D> U[D],O[3],Tmp(&G) ;
00137   for (i=0; i<D; i++) {
00138     U[i].Attach(&G)  ;
00139     U[i].SetRefine() ;
00140   }
00141 
00142   for (i=0; i<3; i++) {
00143     O[i].Attach(&G) ;
00144     O[i].SetBoundaryConditions(BC) ;
00145     O[i].Set(0.0) ;
00146   }
00147 
00148   // Read Velocity fields 
00149   G.ReadSparse(ufile[0]) ;
00150 
00151   for (i=0; i<D; i++) printf("XA/E[%d] %e %e\n",i,G.XA[i],G.XE[i]) ;
00152 
00153   if (!onefile)
00154     for (i=0; i<D; i++) U[i].ReadSparse(ufile[i]) ;
00155   else {
00156     AdaptiveData<D> *UVW[D] ;
00157     for (i=0; i<D; i++) UVW[i]=&U[i] ;
00158 
00159     int which[3]={0,1,2} ;
00160     U[0].ReadSparse(ufile[0],UVW,D,which) ;
00161   }
00162 
00163   for (i=0; i<D; i++) {
00164     U[i].SetBoundaryConditions(BC) ;
00165     printf("U[%d] min/max raw data %e %e\n",i,U[i].Min(), U[i].Max()) ;
00166   }
00167   
00168   // Calculate Omega on Adaptive Grid
00169   int ia, ie ;
00170   if (!absomega) ia=ie=component ;
00171   else { ia= (D==2) ? 2 : 0;  ie=2 ; }   
00172 
00173   for (i=ia ; i<=ie; i++) {
00174     j=(i+1)%3 , k=(i+2)%3 ;
00175     O[i].ApplyOp(BC[j],&U[k] , j , FIRSTDERIVATIVE) ;
00176      Tmp.ApplyOp(BC[k],&U[j] , k , FIRSTDERIVATIVE) ;
00177     O[i].Sub(&O[i] , &Tmp) ;
00178   }
00179  
00180   printf("omega calculated\n") ;
00181 
00182   if (stat) {
00183     assert(D ==3) ;
00184     assert(ia==0) ;
00185     assert(ie==2) ;
00186 
00187     for (i=0; i<3; i++)
00188       for (j=0; j<3; j++ ) O[i].ApplyOp(BC[j],&O[i], j , INVERSETRANSFORM ) ;
00189 
00190     double om[3][3][10000],dm[10000],ens=0,enrg=0 ;
00191     for (i=0; i<3; i++) 
00192       for (j=0; j<3; j++) {
00193         Tmp.Mul(&O[i],&O[j]) ;
00194         for (k=0; k<3; k++) Tmp.ApplyOp(BC[k],&Tmp, k , TRANSFORM ) ;
00195         Tmp.ToUniform(&M,L) ;
00196         for (k=0; k<3; k++) M.ApplyOp(BC[k], &M , k , INVERSETRANSFORM) ;
00197         
00198         M.MeanVariance(om[i][j] , dm , statdir) ;
00199       }
00200 
00201     ens=0 ;
00202     for (i=0; i<=e[1]; i++) {
00203       double r=((double)i)/e[1], y=M.Ext.XE[1]*r+(1-r)*M.Ext.XA[1] ;
00204       printf("stat %e %e %e %e   %e %e %e   %e %e %e\n",y,
00205              om[0][0][i],om[0][1][i],om[0][2][i] , 
00206              om[1][0][i],om[1][1][i],om[1][2][i] , 
00207              om[2][0][i],om[2][1][i],om[2][2][i]) ;
00208 
00209       ens += om[0][0][i]*om[0][0][i] + om[1][1][i]*om[2][2][i] ;
00210     }
00211 
00212     // mean energy and enstrophy
00213     enrg=0 ;
00214     for (i=0; i<3; i++) {
00215       for (j=0; j<3; j++ ) U[i].ApplyOp(BC[j],&U[i], j , INVERSETRANSFORM ) ;
00216       Tmp.Mul(&U[i],&U[i]) ;
00217       for (j=0; j<3; j++ ) Tmp.ApplyOp(BC[j],&Tmp, j , INVERSETRANSFORM ) ;
00218 
00219       enrg += Tmp.Integrate() ;
00220     }
00221 
00222     printf("eng/ens %e %e\n",enrg,ens) ;
00223 
00224     return 0 ;
00225   }
00226 
00227   if (absomega) {
00228 
00229     // Calculate Abs Omega on nonadaptive Grid
00230     for (i=ia; i<=ie; i++) {
00231       O[i].ToUniform(&M,L) ;
00232       for (j=0; j<D; j++) M.ApplyOp(BC[j], &M , j , INVERSETRANSFORM) ;
00233       M.Mul(&M,&M) ;
00234       if (i==ia) W.Copy (&M) ;
00235       else       W.PPlus(&M) ;
00236     }
00237     W.Pow(&W,0.5) ;
00238   }
00239   else {
00240     O[component].ToUniform(&W,L) ;
00241     for (j=0; j<D; j++) W.ApplyOp(BC[j], &W , j , INVERSETRANSFORM) ;
00242   }
00243  
00244   if (smooth) {
00245     for (j=0; j<D; j++) W.ApplyOp(BC[j], &W , j , SMOOTH0) ;
00246   }
00247  
00248   printf("Mix/Max %e %e\n",W.Min(),W.Max()) ;
00249 
00250   if (writeomega) {
00251     if (!restrict) {
00252       W.WriteUDF(ofile,NULL,true) ;
00253     }
00254     else {
00255       // get subset of the full grid
00256       UniformData<D> MM(a,e,&WC) ;
00257       MM.Set(&W) ;
00258       MM.WriteUDF(ofile ,NULL, true) ;
00259     }
00260   }
00261 }
Generated at Mon Aug 19 10:02:32 2002 for AWFD by