analyticalSol.h Source File

00001 #ifndef ANALYTICALSOL_H
00002 #define ANALYTICALSOL_H
00003 
00010 #ifdef SINGLE
00011 #define real float
00012 #else
00013 #define real double
00014 #endif
00015 
00016 extern "C" {
00017 #ifdef SINGLE
00018 void salegpsd_(int *nMax,int *mMax,real *cosTheta,real *pl,real *ps,real *pd,int *ldp);
00019 #else
00020 void alegpsd_(int *nMax,int *mMax,real *cosTheta,real *pl,real *ps,real *pd,int *ldp);
00021 #endif
00022 } 
00023 
00024 
00025 
00026 #include <iostream>
00027 #include <string>
00028 #include <cmath>
00029 #include <vector>
00030 #include "LLSSolver.h"
00031 #include "constant.h"
00032 
00033 using namespace std;
00034 
00042 class analyticalSol {
00043   
00044  public:
00045 
00050   analyticalSol();
00051 
00088   void solve(vector<real> &p,  
00089              const vector<real> &rP,   
00090              const vector<real> &thetaP, 
00091              const vector<real> &phiP,   
00092              vector<real> &omegaX,       
00093              const vector<real> &rOmegaX, 
00094              const vector<real> &thetaOmegaX, 
00095              const vector<real> &phiOmegaX, 
00096              vector<real> &omegaY, 
00097              const vector<real> &rOmegaY, 
00098              const vector<real> &thetaOmegaY, 
00099              const vector<real> &phiOmegaY, 
00100              vector<real> &omegaZ, 
00101              const vector<real> &rOmegaZ, 
00102              const vector<real> &thetaOmegaZ, 
00103              const vector<real> &phiOmegaZ, 
00104              vector<real> &u, 
00105              const vector<real> &rU, 
00106              const vector<real> &thetaU, 
00107              const vector<real> &phiU,
00108              vector<real> &v, 
00109              const vector<real> &rV, 
00110              const vector<real> &thetaV, 
00111              const vector<real> &phiV,
00112              vector<real> &w, 
00113              const vector<real> &rW, 
00114              const vector<real> &thetaW, 
00115              const vector<real> &phiW, 
00116              int nTruncate,  
00117              real ReynoldsNumber,  
00118              vector<real> &solvedABCoeffOneParticle, 
00119              vector<real> &UFixedParticles, 
00120              vector<real> &SFixedParticles, 
00121              vector<real> &VFixedParticles, 
00122              int &fixedParticles, 
00123              int rotatingParticles, 
00124              int particle, 
00125              real **angVel, 
00126              real **angVelAcc, 
00127              real **torqueOld, 
00128              real **torqueNew 
00129              );
00130 
00135   ~analyticalSol();
00136 
00137  private: 
00138 
00139 // modify fields of pressure and vorticity
00140   void modifiedFields(vector<real> &p,
00141                       const vector<real> &rP,
00142                       const vector<real> &thetaP,
00143                       const vector<real> &phiP,
00144                       vector<real> &omegaX,
00145                       const vector<real> &rOmegaX,
00146                       const vector<real> &thetaOmegaX,
00147                       const vector<real> &phiOmegaX,
00148                       vector<real> &omegaY,
00149                       const vector<real> &rOmegaY,
00150                       const vector<real> &thetaOmegaY,
00151                       const vector<real> &phiOmegaY,
00152                       vector<real> &omegaZ,
00153                       const vector<real> &rOmegaZ,
00154                       const vector<real> &thetaOmegaZ,
00155                       const vector<real> &phiOmegaZ,
00156                       int particle, // mth particle
00157                       real **angVel,
00158                       real **angVelAcc);
00159 
00160 // modify fields of velocity due to rotation
00161   void modifiedVelocity(vector<real> &u,
00162                         const vector<real> &rU,
00163                         const vector<real> &thetaU,
00164                         const vector<real> &phiU,
00165                         vector<real> &v,
00166                         const vector<real> &rV,
00167                         const vector<real> &thetaV,
00168                         const vector<real> &phiV,
00169                         vector<real> &w,
00170                         const vector<real> &rW,
00171                         const vector<real> &thetaW,
00172                         const vector<real> &phiW,
00173                         int particle, // mth particle
00174                         real **angVel,
00175                         real **angVelAcc);
00176 
00177   // solving coefficients of A, B, and their tilde counterparts
00178   void solveABCoeff(const vector<real> &p,
00179                     const vector<real> &rP,
00180                     const vector<real> &thetaP,
00181                     const vector<real> &phiP,
00182                     const vector<real> &omegaX,
00183                     const vector<real> &rOmegaX,
00184                     const vector<real> &thetaOmegaX,
00185                     const vector<real> &phiOmegaX,
00186                     const vector<real> &omegaY,
00187                     const vector<real> &rOmegaY,
00188                     const vector<real> &thetaOmegaY,
00189                     const vector<real> &phiOmegaY,
00190                     const vector<real> &omegaZ,
00191                     const vector<real> &rOmegaZ,
00192                     const vector<real> &thetaOmegaZ,
00193                     const vector<real> &phiOmegaZ,
00194                     vector<real> &UFixedParticles,
00195                     vector<real> &SFixedParticles,
00196                     vector<real> &VFixedParticles,
00197                     int &fixedParticles);
00198 
00199   void setB(const vector<real> &p,
00200             const vector<real> &omegaX,
00201             const vector<real> &omegaY,
00202             const vector<real> &omegaZ);
00203 
00204   void setWeightMatrix(const vector<real> &rP,
00205                        const vector<real> &thetaP,
00206                        const vector<real> &phiP,
00207                        const vector<real> &rOmegaX,
00208                        const vector<real> &thetaOmegaX,
00209                        const vector<real> &phiOmegaX,
00210                        const vector<real> &rOmegaY,
00211                        const vector<real> &thetaOmegaY,
00212                        const vector<real> &phiOmegaY,
00213                        const vector<real> &rOmegaZ,
00214                        const vector<real> &thetaOmegaZ,
00215                        const vector<real> &phiOmegaZ);
00216 
00217   void setWeightMatrixFromP(const vector<real> &rP,
00218                             const vector<real> &thetaP,
00219                             const vector<real> &phiP);
00220 
00221   void setWeightMatrixFromOmega(const vector<real> &rOmega,
00222                                 const vector<real> &thetaOmega,
00223                                 const vector<real> &phiOmega,
00224                                 int nStart,
00225                                 int nEnd,
00226                                 int vorticityXYZ);
00227 
00228 
00229   // update velocity
00230   void updateVelocity(vector<real> &u,
00231                       const vector<real> &rU,
00232                       const vector<real> &thetaU,
00233                       const vector<real> &phiU,
00234                       vector<real> &v,
00235                       const vector<real> &rV,
00236                       const vector<real> &thetaV,
00237                       const vector<real> &phiV,
00238                       vector<real> &w,
00239                       const vector<real> &rW,
00240                       const vector<real> &thetaW,
00241                       const vector<real> &phiW);
00242 
00243   void setVelocityWeight(real rVelocity, 
00244                          real thetaVelocity,
00245                          real phiVelocity,
00246                          int velocityXYZ);
00247   
00248   
00249   real sphericalToCartesian(real vectorR,
00250                             real vectorTheta,
00251                             real vectorPhi,
00252                             real theta,
00253                             real phi,
00254                             int vectorXYZ);
00255   
00256   real * allocate1dArray(int size);
00257   void delete1dArray(real * x);
00258 
00259   void allocateMemory();
00260   void releaseMemory(); 
00261 
00262   int i, j, k; // dummy index
00263 
00264   int n;  // degree
00265   int m;  // order
00266 
00267   real rToN;
00268  
00269   // number of pressure points and vortivity points
00270   int nP;
00271   int nOmegaX;
00272   int nOmegaY;
00273   int nOmegaZ;
00274   int nEquations; 
00275 
00276   int Nc;
00277   int nCoeff;
00278 
00279   real Re;
00280 
00281   int nU;
00282   int nV;
00283   int nW;
00284 
00285   real r, theta, phi;
00286 
00287   real cosTheta;
00288   real omegaR, omegaTheta, omegaPhi;
00289 
00290   real * weightMatrix;
00291   real * b;
00292   real * ABCoeff;
00293 
00294   vector<real> velocityWeight;
00295   real uR, uTheta, uPhi;
00296 
00297   // for associated Legendre polynomials
00298   real *pl;
00299   real *ps;
00300   real *pd;
00301   int ldp;
00302 
00303 
00304 };
00305 
00306 #endif