SWE
/import/home/rettenbs/src/SWE/submodules/swe_solvers/src/solver/AugRie.hpp
00001 
00063 /*
00064  * TODO: store nLow/nHigh variables in array[2]
00065  */
00066 
00067 
00068 #ifndef AUGRIE_HPP_
00069 #define AUGRIE_HPP_
00070 
00071 #include "WavePropagation.hpp"
00072 
00073 #include <cassert>
00074 #include <cmath>
00075 #include <iostream>
00076 #include <string>
00077 #include <vector>
00078 
00088 //#define correctrarefactions
00089 //#define complexsteadystatewave
00090 
00091 namespace solver {
00092   template <typename T> class AugRie;
00093 }
00094 
00100 template <typename T> class solver::AugRie: public WavePropagation<T> {
00101   private: //explicit for unit tests
00102     //use nondependent names (template base class)
00103     using solver::WavePropagation<T>::zeroTol;
00104     using solver::WavePropagation<T>::g;
00105     using solver::WavePropagation<T>::dryTol;
00106 
00107     using solver::WavePropagation<T>::hLeft;
00108     using solver::WavePropagation<T>::hRight;
00109     using solver::WavePropagation<T>::huLeft;
00110     using solver::WavePropagation<T>::huRight;
00111     using solver::WavePropagation<T>::bLeft;
00112     using solver::WavePropagation<T>::bRight;
00113     using solver::WavePropagation<T>::uLeft;
00114     using solver::WavePropagation<T>::uRight;
00115 
00116     using solver::WavePropagation<T>::wetDryState;
00117     using solver::WavePropagation<T>::DryDry;
00118     using solver::WavePropagation<T>::WetWet;
00119     using solver::WavePropagation<T>::WetDryInundation;
00120     using solver::WavePropagation<T>::WetDryWall;
00121     using solver::WavePropagation<T>::WetDryWallInundation;
00122     using solver::WavePropagation<T>::DryWetInundation;
00123     using solver::WavePropagation<T>::DryWetWall;
00124     using solver::WavePropagation<T>::DryWetWallInundation;
00125 
00126     using solver::WavePropagation<T>::storeParameters;
00127 
00128 
00130     const T newtonTol;
00132     const int maxNumberOfNewtonIterations;
00133 
00135     T hMiddle;
00136 
00138     T middleStateSpeeds[2];
00139 
00140 
00144     enum RiemannStructure {
00145       DrySingleRarefaction,  
00146       SingleRarefactionDry,  
00147       ShockShock,            
00148       ShockRarefaction,      
00149       RarefactionShock,      
00150       RarefactionRarefaction 
00151     };
00152 
00154     RiemannStructure riemannStructure;
00155 
00156   public:
00166     AugRie(     T i_dryTolerance =                  (T) 0.01,
00167             T i_gravity =                       (T) 9.81,
00168             T i_newtonTolerance =               (T) 0.000001,
00169             int i_maxNumberOfNewtonIterations =     10,
00170             T i_zeroTolerance =                   (T) 0.00001 ):
00171               WavePropagation<T>( i_dryTolerance, i_gravity, i_zeroTolerance ),
00172               newtonTol(i_newtonTolerance),
00173               maxNumberOfNewtonIterations(i_maxNumberOfNewtonIterations) {
00174 
00175 #ifndef NDEBUG
00176 #ifndef SUPPRESS_SOLVER_DEBUG_OUTPUT
00177       //print some information about the used solver
00178       std::cout << "  *** solver::AugRie created" << std::endl
00179                 << "    zeroTolerance=" << zeroTol << std::endl
00180                 << "    gravity=" << g << std::endl
00181                 << "    dryTolerance=" << dryTol << std::endl
00182                 << "    newtonTolerance=" << newtonTol << std::endl
00183                 << "    maxNumberNumberOfNewtonIterations=" << maxNumberOfNewtonIterations << std::endl
00184                 << "\n  optional features:"
00185                 #ifdef correctrarefactions
00186                   << "   correctrarefactions"
00187                 #endif
00188                 #ifdef complexsteadystatewave
00189                   << "   complexsteadystatewave"
00190                 #endif
00191                 << "\n  ***\n\n";
00192 #endif
00193 #endif
00194 
00195     }
00196     //declare variables which are used over and over again
00197 #if 0
00198     T sqrt_g_h[2];
00199     T sqrt_h[2];
00200 
00201 #define sqrt_g_hLeft  (sqrt_g_h[0])
00202 #define sqrt_g_hRight (sqrt_g_h[1])
00203 
00204 #define sqrt_hLeft    (sqrt_h[0])
00205 #define sqrt_hRight   (sqrt_h[1])
00206 
00207 #else
00208     T sqrt_g_hLeft;
00209     T sqrt_g_hRight;
00210 
00211     T sqrt_hLeft;
00212     T sqrt_hRight;
00213 #endif
00214 
00215 
00216     T sqrt_g;
00217 
00218 
00224     void computeNetUpdates ( const T &i_hLeft,  const T &i_hRight,
00225                              const T &i_huLeft, const T &i_huRight,
00226                              const T &i_bLeft,  const T &i_bRight,
00227 
00228                              T &o_hUpdateLeft,
00229                              T &o_hUpdateRight,
00230                              T &o_huUpdateLeft,
00231                              T &o_huUpdateRight,
00232                              T &o_maxWaveSpeed ) {
00233       // store parameters to member variables
00234       storeParameters( i_hLeft, i_hRight,
00235                        i_huLeft, i_huRight,
00236                        i_bLeft, i_bRight );
00237 
00238       //set speeds to zero (will be determined later)
00239       uLeft = uRight = 0.;
00240 
00241       //reset net updates and the maximum wave speed
00242       o_hUpdateLeft = o_hUpdateRight = o_huUpdateLeft = o_huUpdateRight  = (T)0.;
00243       o_maxWaveSpeed = (T)0.;
00244 
00245       //determine the wet/dry state and compute local variables correspondingly
00246       determineWetDryState();
00247 
00248       if (wetDryState == DryDry) //nothing to do in a dry region
00249         return;
00250 
00251       //precompute some terms which are fixed during
00252       //the computation after some specific point
00253       sqrt_g = std::sqrt(g);
00254       sqrt_hLeft = std::sqrt(hLeft);
00255       sqrt_hRight = std::sqrt(hRight);
00256 
00257       sqrt_g_hLeft = sqrt_g*sqrt_hLeft;
00258       sqrt_g_hRight = sqrt_g*sqrt_hRight;
00259 
00260 
00261       //where to store the three waves
00262       T fWaves[3][2];
00263       //and their speeds
00264       T waveSpeeds[3];
00265 
00266       //compute the augmented decomposition
00267       //  (thats the place where the computational work is done..)
00268       computeWaveDecomposition( fWaves,
00269                                 waveSpeeds );
00270 
00271 
00272       //compute the updates from the three propagating waves
00273       //A^-\delta Q = \sum{s[i]<0} \beta[i] * r[i] = A^-\delta Q = \sum{s[i]<0} Z^i
00274       //A^+\delta Q = \sum{s[i]>0} \beta[i] * r[i] = A^-\delta Q = \sum{s[i]<0} Z^i
00275       for (int waveNumber = 0; waveNumber < 3; waveNumber++) {
00276         if (waveSpeeds[waveNumber] < -zeroTol) { //left going
00277           o_hUpdateLeft +=  fWaves[waveNumber][0];
00278           o_huUpdateLeft += fWaves[waveNumber][1];
00279         }
00280 
00281         else if (waveSpeeds[waveNumber] > zeroTol) { //right going
00282           o_hUpdateRight +=  fWaves[waveNumber][0];
00283           o_huUpdateRight += fWaves[waveNumber][1];
00284         }
00285         else { //TODO: this case should not happen mathematically, but it does. Where is the bug? Machine accuracy only?
00286           o_hUpdateLeft +=  (T).5*fWaves[waveNumber][0];
00287           o_huUpdateLeft += (T).5*fWaves[waveNumber][1];
00288 
00289           o_hUpdateRight +=  (T).5*fWaves[waveNumber][0];
00290           o_huUpdateRight += (T).5*fWaves[waveNumber][1];
00291         }
00292         
00293         //no wave speeds => zero strength fWaves
00294 //        assert( std::fabs(fWaves[waveNumber][0]) < zeroTol );
00295 //        assert( std::fabs(fWaves[waveNumber][1]) < zeroTol );
00296       }
00297 
00298 #ifndef NDEBUG
00299       if(wetDryState == DryWetWall)
00300         assert( std::fabs(o_hUpdateLeft) < zeroTol && std::fabs(o_huUpdateLeft) < zeroTol );
00301       else if(wetDryState == WetDryWall)
00302         assert( std::fabs(o_hUpdateRight) < zeroTol && std::fabs(o_huUpdateRight) < zeroTol );
00303 #endif
00304 
00305       //compute maximum wave speed (-> CFL-condition)
00306       waveSpeeds[0] = std::fabs(waveSpeeds[0]);
00307       waveSpeeds[1] = std::fabs(waveSpeeds[1]);
00308       waveSpeeds[2] = std::fabs(waveSpeeds[2]);
00309 
00310       o_maxWaveSpeed = std::max(waveSpeeds[0], waveSpeeds[1]);
00311       o_maxWaveSpeed = std::max(o_maxWaveSpeed, waveSpeeds[2]);
00312     }
00313 
00314   private:
00318     void determineWetDryState() {
00319       //compute speeds or set them to zero (dry cells)
00320       if (hLeft > dryTol) {
00321         uLeft = huLeft / hLeft;
00322       }
00323       else {
00324         bLeft += hLeft;
00325         hLeft = huLeft = uLeft = 0;
00326       }
00327 
00328       if (hRight > dryTol) {
00329         uRight = huRight / hRight;
00330       }
00331       else {
00332         bRight += hRight;
00333         hRight = huRight = uRight = 0;
00334       }
00335 
00336       //test for simple wet/wet case since this is most probably the
00337       //most frequently executed branch
00338       if (hLeft >= dryTol && hRight >= dryTol) {
00339         wetDryState =  WetWet;
00340       }
00341 
00342       //check for the dry/dry-case
00343       else if (hLeft < dryTol && hRight < dryTol)
00344         wetDryState =  DryDry;
00345 
00346       //we have a shoreline: one cell dry, one cell wet
00347 
00348       //check for simple inundation problems
00349       // (=> dry cell lies lower than the wet cell)
00350       else if (hLeft < dryTol && hRight + bRight > bLeft)
00351         wetDryState =  DryWetInundation;
00352 
00353       else if (hRight < dryTol && hLeft + bLeft > bRight)
00354         wetDryState =  WetDryInundation;
00355 
00356       //dry cell lies higher than the wet cell
00357       else if (hLeft < dryTol) {
00358         //lets check if the momentum is able to overcome the difference in height
00359         //  => solve homogeneous Riemann-problem to determine the middle state height
00360         //     which would arise if there is a wall (wall-boundary-condition)
00361         //       \cite[ch. 6.8.2]{george2006finite})
00362         //       \cite[ch. 5.2]{george2008augmented}
00363         computeMiddleState( hRight,  hRight,
00364                            -uRight, uRight,
00365                            -huRight, huRight,
00366                            maxNumberOfNewtonIterations );
00367 
00368         if (hMiddle + bRight > bLeft) {
00369           //momentum is large enough, continue with the original values
00370 //          bLeft = hMiddle + bRight;
00371           wetDryState =  DryWetWallInundation;
00372         }
00373         else {
00374           //momentum is not large enough, use wall-boundary-values
00375           hLeft = hRight;
00376           uLeft = -uRight;
00377           huLeft = -huRight;
00378           bLeft = bRight = (T)0.;
00379           wetDryState =  DryWetWall;
00380         }
00381       }
00382 
00383       else if (hRight < dryTol) {
00384         //lets check if the momentum is able to overcome the difference in height
00385         //  => solve homogeneous Riemann-problem to determine the middle state height
00386         //     which would arise if there is a wall (wall-boundary-condition)
00387         //       \cite[ch. 6.8.2]{george2006finite})
00388         //       \cite[ch. 5.2]{george2008augmented}
00389         computeMiddleState( hLeft,  hLeft,
00390                             uLeft, -uLeft,
00391                             huLeft, -huLeft,
00392                             maxNumberOfNewtonIterations );
00393 
00394         if (hMiddle + bLeft > bRight) {
00395           //momentum is large enough, continue with the original values
00396 //          bRight = hMiddle + bLeft;
00397           wetDryState =  WetDryWallInundation;
00398         }
00399         else {
00400           hRight = hLeft;
00401           uRight = -uLeft;
00402           huRight = -huLeft;
00403           bRight = bLeft = (T)0.;
00404           wetDryState =  WetDryWall;
00405         }
00406       }
00407       //done with all cases
00408       else assert(false);
00409 
00410       //limit the effect of the source term if there is a "wall"
00411       //\cite[end of ch. 5.2?]{george2008augmented}
00412       //\cite[rpn2ez_fast_geo.f][levequeclawpack]
00413       if ( wetDryState == DryWetWallInundation )
00414         bLeft = hRight + bRight;
00415 
00416       else if( wetDryState == WetDryWallInundation )
00417         bRight = hLeft + bLeft;
00418     }
00419 
00426     inline void computeWaveDecomposition( T o_fWaves[3][2],
00427                                           T o_waveSpeeds[3] )   {
00428       //compute eigenvalues of the jacobian matrices in states Q_{i-1} and Q_{i} (char. speeds)
00429       T characteristicSpeeds[2];
00430       characteristicSpeeds[0] = uLeft - sqrt_g_hLeft;
00431       characteristicSpeeds[1] = uRight + sqrt_g_hRight;
00432 
00433       //compute "Roe speeds"
00434       T hRoe = (T).5 * (hRight + hLeft);
00435       T uRoe = uLeft * sqrt_hLeft + uRight * sqrt_hRight;
00436       uRoe /= sqrt_hLeft + sqrt_hRight;
00437 
00438       T roeSpeeds[2];
00439       //optimization for dumb compilers
00440       T sqrt_g_hRoe = std::sqrt(g*hRoe);
00441       roeSpeeds[0] = uRoe - sqrt_g_hRoe;
00442       roeSpeeds[1] = uRoe + sqrt_g_hRoe;
00443 
00444       //compute the middle state of the homogeneous Riemann-Problem
00445       if(wetDryState != WetDryWall && wetDryState != DryWetWall) { //case WDW and DWW was computed in determineWetDryState already
00446         computeMiddleState( hLeft, hRight,
00447                             uLeft, uRight,
00448                             huLeft, huRight
00449         );
00450       }
00451 
00452       //compute extended eindfeldt speeds (einfeldt speeds + middle state speeds)
00453       //  \cite[ch. 5.2]{george2008augmented}, \cite[ch. 6.8]{george2006finite}
00454       T extEinfeldtSpeeds[2]= {(T)0, (T)0};
00455       if( wetDryState == WetWet ||
00456           wetDryState == WetDryWall ||
00457           wetDryState == DryWetWall ) {
00458         extEinfeldtSpeeds[0] = std::min(characteristicSpeeds[0], roeSpeeds[0]);
00459         extEinfeldtSpeeds[0] = std::min(extEinfeldtSpeeds[0], middleStateSpeeds[1]);
00460 
00461         extEinfeldtSpeeds[1] = std::max(characteristicSpeeds[1], roeSpeeds[1]);
00462         extEinfeldtSpeeds[1] = std::max(extEinfeldtSpeeds[1], middleStateSpeeds[0]);
00463       }
00464       else if (hLeft < dryTol) { //ignore undefined speeds
00465         extEinfeldtSpeeds[0] = std::min(roeSpeeds[0], middleStateSpeeds[1]);
00466         extEinfeldtSpeeds[1] = std::max(characteristicSpeeds[1], roeSpeeds[1]);
00467 
00468         assert(middleStateSpeeds[0] < extEinfeldtSpeeds[1]);
00469       }
00470       else if (hRight < dryTol) { //ignore undefined speeds
00471         extEinfeldtSpeeds[0] = std::min(characteristicSpeeds[0], roeSpeeds[0]);
00472         extEinfeldtSpeeds[1] = std::max(roeSpeeds[1], middleStateSpeeds[0]);
00473 
00474         assert(middleStateSpeeds[1] > extEinfeldtSpeeds[0]);
00475       }
00476       else {
00477         assert(false);
00478       }
00479 
00480       //HLL middle state
00481       //  \cite[theorem 3.1]{george2006finite}, \cite[ch. 4.1]{george2008augmented}
00482       T hLLMiddleHeight = huLeft - huRight + extEinfeldtSpeeds[1] * hRight - extEinfeldtSpeeds[0] * hLeft;
00483       hLLMiddleHeight /= extEinfeldtSpeeds[1] - extEinfeldtSpeeds[0];
00484       hLLMiddleHeight = std::max(hLLMiddleHeight, (T)0.);
00485 
00486       //define eigenvalues
00487       T eigenValues[3];
00488       eigenValues[0] = extEinfeldtSpeeds[0];
00489       //eigenValues[1] --> corrector wave
00490       eigenValues[2] = extEinfeldtSpeeds[1];
00491 
00492       //define eigenvectors
00493       T eigenVectors[3][3];
00494 
00495       //set first and third eigenvector
00496       eigenVectors[0][0] = (T)1;
00497       eigenVectors[0][2] = (T)1;
00498 
00499       eigenVectors[1][0] = eigenValues[0];
00500       eigenVectors[1][2] = eigenValues[2];
00501 
00502       eigenVectors[2][0] = eigenValues[0]*eigenValues[0];
00503       eigenVectors[2][2] = eigenValues[2]*eigenValues[2];
00504 
00505 
00506       //compute rarefaction corrector wave
00507       //  \cite[ch. 6.7.2]{george2006finite}, \cite[ch. 5.1]{george2008augmented}
00508       bool strongRarefaction = false;
00509 #ifdef correctrarefactions
00510       if ( (riemannStructure == ShockRarefaction ||
00511            riemannStructure == RarefactionShock ||
00512            riemannStructure == RarefactionRarefaction) &&
00513            hMiddle > dryTol ) { //limit to ensure non-negative depth
00514 
00515         //TODO: GeoClaw, riemann_aug_JCP; No explicit boundaries for "strong rarefaction" in literature?
00516         T rareBarrier[2] = {.5, .9};
00517 
00518         //lower rare barrier in the case of a transsonic rarefaction
00519         if ( (riemannStructure == RarefactionShock || riemannStructure == RarefactionRarefaction) &&
00520              eigenValues[0]*middleStateSpeeds[0] < (T).0 ) { //transsonic rarefaction, first wave family
00521           rareBarrier[0] = (T).2;
00522         }
00523         else if ( (riemannStructure == ShockRarefaction || riemannStructure == RarefactionRarefaction) &&
00524                   eigenValues[2]*middleStateSpeeds[1] < (T)0. ) { //transsonic rarefaction, second wave family
00525           rareBarrier[0] = (T).2;
00526         }
00527 
00528 
00529         T sqrt_g_hMiddle = std::sqrt(g*hMiddle);
00530         //determine the max. rarefaction size (distance between head and tail)
00531         T rareFactionSize[2];
00532         rareFactionSize[0] = (T)3. * (sqrt_g_hLeft - sqrt_g_hMiddle);
00533         rareFactionSize[1] = (T)3. * (sqrt_g_hRight - sqrt_g_hMiddle);
00534 
00535         T maxRareFactionSize = std::max(rareFactionSize[0], rareFactionSize[1]);
00536 
00537         //set the eigenvalue of the corrector wave in the case of a "strong rarefaction"
00538         if ( maxRareFactionSize > rareBarrier[0] * (eigenValues[2] - eigenValues[0]) &&
00539              maxRareFactionSize < rareBarrier[1] * (eigenValues[2] - eigenValues[0])
00540         ) {
00541           strongRarefaction = true;
00542           if (rareFactionSize[0] > rareFactionSize[1])
00543             eigenValues[1] = middleStateSpeeds[0];
00544           else
00545             eigenValues[1] = middleStateSpeeds[1];
00546         }
00547 
00548         //TODO: implemented in clawpack, why?
00549 //        if ( hMiddle < std::min(hLeft, hRight)/5.) { //middle state in an HLL solve
00550 //          strongRarefaction = false;
00551 //        }
00552       }
00553 #endif
00554 
00555       //set 2nd eigenvector
00556       if(strongRarefaction == false) {
00557         eigenValues[1] = (T).5*(eigenValues[0]+eigenValues[2]);
00558         eigenVectors[0][1] = 0.;
00559         eigenVectors[1][1] = 0.;
00560         eigenVectors[2][1] = 1.;
00561       }
00562       else {
00563         eigenVectors[0][1] = (T)1.;
00564         eigenVectors[1][1] = eigenValues[1];
00565         eigenVectors[2][1] = eigenValues[1]*eigenValues[1];
00566       }
00567 
00568       //compute the jump in state
00569       T rightHandSide[3];
00570       rightHandSide[0] = hRight - hLeft;
00571       rightHandSide[1] = huRight - huLeft;
00572       rightHandSide[2] = huRight * uRight + (T).5 * g * hRight * hRight
00573                       -( huLeft  * uLeft  + (T).5 * g * hLeft  * hLeft );
00574 
00575       //compute steady state wave
00576       //  \cite[ch. 4.2.1 \& app. A]{george2008augmented}, \cite[ch. 6.2 \& ch. 4.4]{george2006finite}
00577       T steadyStateWave[2];
00578       T hBar = (hLeft + hRight)*(T).5;
00579 #ifdef complexsteadystatewave
00580       T lLambdaBar = (uLeft + uRight)*(uLeft + uRight)*.25 - g*hBar;
00581 
00582       T lLambdaTilde = std::max((T)0., uLeft*uRight) - g*hBar;
00583 
00584       //near sonic as defined in geoclaw (TODO: literature?)
00585       if( ( std::fabs(lLambdaBar) < zeroTol )                                             ||
00586           ( lLambdaBar*lLambdaTilde < zeroTol )                                           ||
00587           ( lLambdaBar*eigenValues[0]*eigenValues[1] < zeroTol )                          ||
00588           ( std::min( std::fabs( eigenValues[0]), std::fabs(eigenValues[2]) ) < zeroTol ) ||
00589           ( eigenValues[0] < (T)0. && middleStateSpeeds[0] > (T)0. )                      ||
00590           ( eigenValues[2] > (T)0. && middleStateSpeeds[1] < (T)0. )                      ||
00591           ( (uLeft + sqrt(g*hLeft)) * (uRight+sqrt(g*hRight)) < (T)0. )                   ||
00592           ( (uLeft - sqrt(g*hLeft)) * (uRight-sqrt(g*hRight)) < (T)0. ) ) {
00593 #endif
00594         steadyStateWave[0] = -(bRight - bLeft);
00595         steadyStateWave[1] = -g * hBar * (bRight - bLeft);
00596 #ifdef complexsteadystatewave
00597       }
00598       else {
00599         steadyStateWave[0] = g * (hBar / lLambdaBar) * (bRight - bLeft);
00600         T hTilde = hBar * lLambdaTilde / lLambdaBar;
00601 
00602         //set bounds for problems far from steady state
00603         hTilde = std::max(std::min(hLeft, hRight), hTilde);
00604         hTilde = std::min(std::max(hLeft, hRight), hTilde);
00605         steadyStateWave[1] = -g*hTilde*(bRight - bLeft);
00606       }
00607 #endif
00608 
00609       //preserve depth-positivity
00610       //  \cite[ch. 6.5.2]{george2006finite}, \cite[ch. 4.2.3]{george2008augmented}
00611       if(eigenValues[0] < -zeroTol && eigenValues[2] > zeroTol) { //subsonic
00612         steadyStateWave[0] = std::max(steadyStateWave[0], hLLMiddleHeight*(eigenValues[2]-eigenValues[0])/eigenValues[0]);
00613         steadyStateWave[0] = std::min(steadyStateWave[0], hLLMiddleHeight*(eigenValues[2]-eigenValues[0])/eigenValues[2]);
00614       }
00615       else if(eigenValues[0] > zeroTol) { //supersonic right TODO: motivation?
00616         steadyStateWave[0] = std::max(steadyStateWave[0], -hLeft);
00617         steadyStateWave[0] = std::min(steadyStateWave[0], hLLMiddleHeight*(eigenValues[2]-eigenValues[0])/eigenValues[0]);
00618       }
00619       else if(eigenValues[2] < -zeroTol) { //supersonic left TODO: motivation?
00620         steadyStateWave[0] = std::max(steadyStateWave[0], hLLMiddleHeight*(eigenValues[2]-eigenValues[0])/eigenValues[2]);
00621         steadyStateWave[0] = std::min(steadyStateWave[0], hRight);
00622       }
00623 
00624       //Limit the effect of the source term
00625       //  \cite[ch. 6.4.2]{george2006finite}
00626       steadyStateWave[1] = std::min(steadyStateWave[1], g*std::max(-hLeft*(bRight-bLeft) , -hRight*(bRight-bLeft)));
00627       steadyStateWave[1] = std::max(steadyStateWave[1], g*std::min(-hLeft*(bRight-bLeft) , -hRight*(bRight-bLeft)));
00628 
00629       //no source term in the case of a wall
00630       if( wetDryState == WetDryWall || wetDryState == DryWetWall ) {
00631         assert(std::fabs(steadyStateWave[0])<zeroTol);
00632         assert(std::fabs(steadyStateWave[1])<zeroTol);
00633       }
00634 
00635       rightHandSide[0] -= steadyStateWave[0];
00636       //rightHandSide[1]: no source term
00637       rightHandSide[2] -= steadyStateWave[1];
00638 
00639       //everything is ready, solve the equations!
00640       T beta[3];
00641       solveLinearEquation(eigenVectors, rightHandSide, beta);
00642 
00643       //compute f-waves and wave-speeds
00644       if (wetDryState == WetDryWall) { //zero ghost updates (wall boundary)
00645         //care about the left going wave (0) only
00646         o_fWaves[0][0] = beta[0] * eigenVectors[1][0];
00647         o_fWaves[0][1] = beta[0] * eigenVectors[2][0];
00648 
00649         //set the rest to zero
00650         o_fWaves[1][0] = o_fWaves[1][1] = (T)0.;
00651         o_fWaves[2][0] = o_fWaves[2][1] = (T)0.;
00652 
00653         o_waveSpeeds[0] = eigenValues[0];
00654         o_waveSpeeds[1] = o_waveSpeeds[2] = (T)0.;
00655 
00656         assert(eigenValues[0] < zeroTol);
00657       }
00658       else if (wetDryState == DryWetWall) { //zero ghost updates (wall boundary)
00659         //care about the right going wave (2) only
00660         o_fWaves[2][0] = beta[2] * eigenVectors[1][2];
00661         o_fWaves[2][1] = beta[2] * eigenVectors[2][2];
00662 
00663         //set the rest to zero
00664         o_fWaves[0][0] = o_fWaves[0][1] = (T)0.;
00665         o_fWaves[1][0] = o_fWaves[1][1] = (T)0.;
00666 
00667         o_waveSpeeds[2] = eigenValues[2];
00668         o_waveSpeeds[0] = o_waveSpeeds[1] = 0.;
00669 
00670         assert(eigenValues[2] > -zeroTol);
00671       }
00672       else {
00673         //compute f-waves (default)
00674         for (int waveNumber = 0; waveNumber < 3; waveNumber++) {
00675           o_fWaves[waveNumber][0] = beta[waveNumber] * eigenVectors[1][waveNumber];  //select 2nd and
00676           o_fWaves[waveNumber][1] = beta[waveNumber] * eigenVectors[2][waveNumber];  //3rd component of the augmented decomposition
00677         }
00678 
00679         o_waveSpeeds[0] = eigenValues[0];
00680         o_waveSpeeds[1] = eigenValues[1];
00681         o_waveSpeeds[2] = eigenValues[2];
00682       }
00683     }
00684 
00697     inline void computeMiddleState( const T &i_hLeft, const T &i_hRight,
00698                                     const T &i_uLeft, const T &i_uRight,
00699                                     const T &i_huLeft, const T &i_huRight,
00700                                     const int &i_maxNumberOfNewtonIterations = 1 ) {
00701 
00702       //set everything to zero
00703       hMiddle = (T)0.;
00704       middleStateSpeeds[0] = (T)0.;
00705       middleStateSpeeds[1] = (T)0.;
00706 
00707       //compute local square roots
00708       //(not necessarily the same ones as in computeNetUpdates!)
00709       T l_sqrt_g_hRight = std::sqrt(g*i_hRight);
00710       T l_sqrt_g_hLeft = std::sqrt(g*i_hLeft);
00711 
00712       //single rarefaction in the case of a wet/dry interface
00713       if (i_hLeft < dryTol) {
00714         middleStateSpeeds[1] = middleStateSpeeds[0] = i_uRight - (T)2*l_sqrt_g_hRight;
00715         riemannStructure = DrySingleRarefaction;
00716         return;
00717       }
00718       else if (i_hRight < dryTol) {
00719         middleStateSpeeds[0] = middleStateSpeeds[1] = i_uLeft + (T)2*l_sqrt_g_hLeft;
00720         riemannStructure = SingleRarefactionDry;
00721         return;
00722       }
00723 
00724       //determine the wave structure of the Riemann-problem
00725       riemannStructure = determineRiemannStructure( i_hLeft, i_hRight,
00726                                                     i_uLeft, i_uRight );
00727 
00728       //will be computed later
00729       T sqrt_g_hMiddle = 0.;
00730 
00731       if (riemannStructure == ShockShock) {
00732         /*compute All-Shock Riemann Solution
00733          *  \cite[ch. 13.7]{leveque2002finite}
00734          *
00735          *compute middle state h_m
00736          * => Solve non-linear scalar equation by Newton's method
00737          *    u_l - (h_m - h_l) \sqrt{\frac{g}{2} \left( \frac{1}{h_m} + \frac{1}{h_l} \right) }
00738          * =  u_r + (h_m - h_r) \sqrt{\frac{g}{2} \left( \frac{1}{h_m} + \frac{1}{h_r} \right) }
00739          *
00740          * Therefore determine the root of phi_{ss}:
00741          *\begin{equation}
00742          *  \begin{matrix}
00743          *    \phi_{ss}(h) =&&  u_r - u_l +
00744          *                  (h-h_l) \sqrt{
00745          *                            \frac{g}{2} \left( \frac{1}{h} + \frac{1}{h_l} \right)
00746          *                          } \\
00747          *             &&  +(h-h_r) \sqrt{
00748          *                            \frac{g}{2} \left( \frac{1}{h} + \frac{1}{h_r} \right)
00749          *                          }
00750          *  \end{matrix}
00751          *\end{equation}
00752          *
00753          *\begin{equation}
00754          *  \begin{matrix}
00755          *    \phi_{ss}'(h) = &&\sqrt{\frac{g}{2} \left( \frac{1}{h} + \frac{1}{h_l} \right) }
00756          *                      +\sqrt{\frac{g}{2} \left( \frac{1}{h} + \frac{1}{h_r} \right) }\\
00757          *                    && -\frac{g}{4}
00758          *                        \frac{h-h_l}
00759          *                        {
00760          *                         h^2\sqrt{\frac{g}{2} \left( \frac{1}{h} + \frac{1}{h_l} \right) }
00761          *                        }-
00762          *                        \frac{g}{4}
00763          *                        \frac{h-h_r}
00764          *                        {
00765          *                         h^2\sqrt{\frac{g}{2} \left( \frac{1}{h} + \frac{1}{h_r} \right) }
00766          *                        }
00767          *  \end{matrix}
00768          *\end{equation}
00769          */
00770 
00771 //        hMiddle = (hLeft + hRight)*(T)0.618; //first estimate
00772         hMiddle = std::min(i_hLeft, i_hRight);
00773 
00774         T l_sqrtTermH[2] = {0,0};
00775 
00776         for (int i = 0; i < i_maxNumberOfNewtonIterations; i++) {
00777 //          sqrtTermHLow = std::sqrt((T).5 * g * ((T)1/hMiddle + (T)1/hLeft));
00778           l_sqrtTermH[0] = std::sqrt((T).5 * g * ((hMiddle+i_hLeft)/(hMiddle*i_hLeft)));
00779 //          sqrtTermHHigh = std::sqrt((T).5 * g * ((T)1/hMiddle + (T)1/hRight));
00780           l_sqrtTermH[1] = std::sqrt((T).5 * g * ((hMiddle+i_hRight)/(hMiddle*i_hRight)));
00781 
00782           T phi = i_uRight - i_uLeft + (hMiddle - i_hLeft) * l_sqrtTermH[0] + (hMiddle - i_hRight) * l_sqrtTermH[1];
00783 
00784           if (std::fabs(phi) < newtonTol)
00785             break;
00786 
00787           T derivativePhi = l_sqrtTermH[0] + l_sqrtTermH[1]
00788                           - (T).25*g*(hMiddle - i_hLeft)  / (l_sqrtTermH[0]*hMiddle*hMiddle)
00789                           - (T).25*g*(hMiddle - i_hRight) / (l_sqrtTermH[1]*hMiddle*hMiddle);
00790 
00791           hMiddle = hMiddle - phi/derivativePhi; //Newton step
00792           assert(hMiddle >= dryTol);
00793 
00794 //          if(i == i_maxNumberOfNewtonIterations-1) {
00795 //            std::cerr << "Newton-Method did not converge" << std::endl;
00796 //            std::cerr << "std::fabs(phi): " << std::fabs(phi) << std::endl;
00797 //            std::cerr << "hMiddle: " << hMiddle << std::endl;
00798 //            assert(false);
00799 //          }
00800         }
00801 
00802         sqrt_g_hMiddle = std::sqrt(g*hMiddle);
00803 
00804         //compute middle speed u_m
00805         //\begin{equation}
00806         //  \label{eq:hmshock1stfamsimp}
00807         //  u_m = u_l - (h_m - h_l) \sqrt{\frac{g}{2} \left( \frac{1}{h_m} + \frac{1}{h_l} \right) }
00808         //\end{equation}
00809         //
00810         //\begin{equation}
00811         //  \label{eq:hmshock2ndfamsimp}
00812         //  u_m = u_r + (h_m - h_r) \sqrt{\frac{g}{2} \left( \frac{1}{h_m} + \frac{1}{h_r} \right) }
00813         //\end{equation}
00814 
00815 //        T uMiddleEstimates[2];
00816 //        uMiddleEstimates[0] = uLeft - (hMiddle - hLeft) * l_sqrtTermH[0];
00817 //        uMiddleEstimates[1] = i_uRight + (hMiddle - hRight) * l_sqrtTermH[1];
00818 //        uMiddle = (T).5 * (uMiddleEstimates[0] + uMiddleEstimates[1]);
00819 
00820         //middle state speeds as they are implemented in clawpack, TODO: why?
00821 //        middleStateSpeeds[0] = uMiddleEstimates[0] - sqrt_g_hMiddle;
00822 //        middleStateSpeeds[1] = uMiddleEstimates[1] + sqrt_g_hMiddle;
00823       }
00824 
00825       if (riemannStructure == RarefactionRarefaction) {
00826         //Compute All-Rarefaction Riemann Solution
00827         //  \cite[ch. 13.8.6]{leveque2002finite}
00828 
00829         //compute middle state height h_m
00830         hMiddle = std::max((T)0., i_uLeft - i_uRight + (T)2.*(l_sqrt_g_hLeft + l_sqrt_g_hRight)); //std::max -> Text after formula (13.56), page 279
00831         hMiddle = (T)1/((T)16*g)* hMiddle * hMiddle;
00832 
00833         sqrt_g_hMiddle = std::sqrt(g*hMiddle);
00834 
00835         //middle state speeds as they are implemented in clawpack, why?
00836 //        middleStateSpeeds[0] = uLeft + (T)2.*l_sqrt_g_hLow - (T)3.*sqrt_g_hMiddle;
00837 //        middleStateSpeeds[1] = i_uRight - (T)2.*l_sqrt_g_hHigh + (T)3.*sqrt_g_hMiddle;
00838       }
00839 
00840       if (riemannStructure == ShockRarefaction || riemannStructure == RarefactionShock) {
00841         //compute dam-break Riemann-solution
00842         //TODO: reference
00843         T hMin, hMax;
00844         if (riemannStructure == ShockRarefaction) {
00845           hMin = i_hLeft;
00846           hMax = i_hRight;
00847         }
00848         else {
00849           hMin = i_hRight;
00850           hMax = i_hLeft;
00851         }
00852 
00853         //hMiddle = (hLeft + hRight)*(T)0.618; //first estimate
00854         hMiddle = hMin;
00855 
00856         sqrt_g_hMiddle = std::sqrt(g*hMiddle);
00857         T sqrt_g_hMax = std::sqrt(g*hMax);
00858         for (int i = 0; i < i_maxNumberOfNewtonIterations; i++) {
00859           /*
00860            *compute middle state h_m
00861            * => Solve non-linear scalar equation by Newton's method
00862            *
00863            * Therefore determine the root of phi_{sr}/phi_{rs}:
00864            *\begin{equation}
00865            *  \begin{matrix}
00866            *    h_{min} = \min(h_l, h_r)\\
00867            *    h_{max} = \max(h_l, h_r)
00868            *  \end{matrix}
00869            *\end{equation}
00870            *\begin{equation}
00871            *  \begin{matrix}
00872            *    \phi_{sr}(h) = \phi_{rs}(h) =&& u_r - u_l + (h - h_{min}) \sqrt{
00873            *                                            \frac{g}{2} \left( \frac{1}{h} + \frac{1}{h_{min}} \right)
00874            *                                          } \\
00875            *               && + 2 (\sqrt{gh} - \sqrt{gh_{max}})
00876            *  \end{matrix}
00877            *\end{equation}
00878            *\begin{equation}
00879            *    \phi'_{sr}(h) = \phi'_{rs}(h) = \sqrt{\frac{g}{2} \left( \frac{1}{h} + \frac{1}{h_{min}} \right)  }
00880            *                     -\frac{g}{4} \frac{h-h_{min}}
00881            *                        { h^2
00882            *                          \sqrt{
00883            *                            \frac{g}{2}  \left(\frac{1}{h} + \frac{1}{h_{min}} \right)
00884            *                          }
00885            *                        }
00886            *                     + \sqrt{\frac{g}{h}}
00887            *\end{equation}
00888            */
00889 
00890           T sqrtTermHMin = std::sqrt( (T).5*g* ((hMiddle+hMin)/(hMiddle*hMin)) );
00891 
00892           T phi = i_uRight - i_uLeft + (hMiddle - hMin) * sqrtTermHMin + (T)2*(sqrt_g_hMiddle - sqrt_g_hMax);
00893 
00894           if (std::fabs(phi) < newtonTol)
00895             break;
00896 
00897           T derivativePhi = sqrtTermHMin - (T).25*g * (hMiddle - hMin) / (hMiddle * hMiddle * sqrtTermHMin) + sqrt_g/sqrt_g_hMiddle;
00898 
00899           hMiddle = hMiddle - phi/derivativePhi; //Newton step
00900 
00901           #ifndef NDEBUG
00902           if (hMiddle < hMin) {
00903             std::cout << phi << std::endl;
00904             std::cout << derivativePhi << std::endl;
00905             std::cerr << "hMiddle(" << hMiddle << ") < hMin(" << hMin << ")" << std::endl;
00906             assert(false);
00907           }
00908           #endif
00909 
00910 //          if(i == i_maxNumberOfNewtonIterations-1) {
00911 //            std::cerr << "Newton-Method did not converge" << std::endl;
00912 //            std::cerr << "std::fabs(phi): " << std::fabs(phi) << std::endl;
00913 //            std::cerr << "hMiddle: " << hMiddle << std::endl;
00914 //            assert(false);
00915 //          }
00916 
00917           sqrt_g_hMiddle = std::sqrt(g*hMiddle);
00918         }
00919 
00920 //        //middle state speeds as they are implemented in clawpack, TODO: why?
00921 //        if (riemannStructure == ShockRarefaction) {
00922 //          uMiddle = i_uRight - (T)2.*( l_sqrt_g_hHigh - sqrt_g_hMiddle );
00923 //          middleStateSpeeds[0] = i_uRight - (T)2.*l_sqrt_g_hHigh + sqrt_g_hMiddle;
00924 //          middleStateSpeeds[1] = i_uRight - (T)2.*l_sqrt_g_hHigh + (T)3.*sqrt_g_hMiddle;
00925 //        }
00926 //        else {
00927 //          uMiddle = uLeft + (T)2.*( l_sqrt_g_hLow - sqrt_g_hMiddle );
00928 //          middleStateSpeeds[0] = uLeft + (T)2.*l_sqrt_g_hLow - (T)3.*sqrt_g_hMiddle;
00929 //          middleStateSpeeds[1] = uLeft + (T)2.*l_sqrt_g_hLow - sqrt_g_hMiddle;
00930 //        }
00931       }
00932 
00933       middleStateSpeeds[0] = i_uLeft + (T)2.*l_sqrt_g_hLeft - (T)3.*sqrt_g_hMiddle;
00934       middleStateSpeeds[1] = i_uRight - (T)2.*l_sqrt_g_hRight + (T)3.*sqrt_g_hMiddle;
00935 
00936       assert(hMiddle >= 0);
00937     }
00938 
00950     inline RiemannStructure determineRiemannStructure( const T &i_hLeft, const T &i_hRight,
00951                                                        const T &i_uLeft, const T &i_uRight ) const {
00952       T hMin = std::min(i_hLeft, i_hRight);
00953       T hMax = std::max(i_hLeft, i_hRight);
00954 
00955       T uDif = i_uRight - i_uLeft;
00956 
00957       if( 0 <= (T)2.0*(std::sqrt(g*hMin) - std::sqrt(g*hMax)) + uDif)
00958         return RarefactionRarefaction;
00959 
00960       if( (hMax - hMin) * std::sqrt(g*(T)0.5*( 1/hMax + 1/hMin )) + uDif <= 0 )
00961         return ShockShock;
00962 
00963       if( i_hLeft < i_hRight )
00964         return ShockRarefaction;
00965 
00966       return RarefactionShock;
00967     }
00968 
00977     static inline void solveLinearEquation( const T i_matrix[3][3],
00978                                             const T i_b[3],
00979                                             T o_x[3] ) {
00980 //#if 1
00981       // compute inverse of 3x3 matrix
00982       const T m[3][3] =
00983         {       {        (i_matrix[1][1]*i_matrix[2][2] - i_matrix[1][2]*i_matrix[2][1]), -(i_matrix[0][1]*i_matrix[2][2] - i_matrix[0][2]*i_matrix[2][1]),  (i_matrix[0][1]*i_matrix[1][2] - i_matrix[0][2]*i_matrix[1][1]) },
00984         {   -(i_matrix[1][0]*i_matrix[2][2] - i_matrix[1][2]*i_matrix[2][0]),  (i_matrix[0][0]*i_matrix[2][2] - i_matrix[0][2]*i_matrix[2][0]), -(i_matrix[0][0]*i_matrix[1][2] - i_matrix[0][2]*i_matrix[1][0]) },
00985         {    (i_matrix[1][0]*i_matrix[2][1] - i_matrix[1][1]*i_matrix[2][0]), -(i_matrix[0][0]*i_matrix[2][1] - i_matrix[0][1]*i_matrix[2][0]),  (i_matrix[0][0]*i_matrix[1][1] - i_matrix[0][1]*i_matrix[1][0]) }
00986       };
00987       T d = (i_matrix[0][0]*m[0][0] + i_matrix[0][1]*m[1][0] + i_matrix[0][2]*m[2][0]);
00988 
00989 #ifndef NDEBUG
00990       if (std::fabs(d) < 0.000000001) {
00991         std::cerr << "Division close to zero!" << std::endl;
00992         std::cout << "Matrix: " << std::endl;
00993         std::cout << i_matrix[0][0] << ", " << i_matrix[0][1] << ", " << i_matrix[0][2] << std::endl;
00994         std::cout << i_matrix[1][0] << ", " << i_matrix[1][1] << ", " << i_matrix[1][2] << std::endl;
00995         std::cout << i_matrix[2][0] << ", " << i_matrix[2][1] << ", " << i_matrix[2][2] << std::endl;
00996         std::cout << "b: " << std::endl;
00997         std::cout << i_b[0] << ", " << i_b[1] << ", " << i_b[2] << std::endl;
00998         std::cout << "d: " << d << std::endl;
00999         assert(false);
01000       }
01001 #endif
01002 
01003       // m stores not really the inverse matrix, but the inverse multiplied by d
01004       T s = (T)1/d;
01005 
01006       // compute m*i_b
01007       o_x[0] = (m[0][0]*i_b[0] + m[0][1]*i_b[1] + m[0][2]*i_b[2])*s;
01008       o_x[1] = (m[1][0]*i_b[0] + m[1][1]*i_b[1] + m[1][2]*i_b[2])*s;
01009       o_x[2] = (m[2][0]*i_b[0] + m[2][1]*i_b[1] + m[2][2]*i_b[2])*s;
01010 
01011 //      return;
01012 //#else
01013 //      T origDet = computeDeterminant(i_matrix);
01014 //
01015 //      if (std::fabs(origDet) > zeroTol) {
01016 //        T modifiedMatrix[3][3];
01017 //
01018 //        for (int column = 0; column < 3; column++)
01019 //        {
01020 //          memcpy(modifiedMatrix, i_matrix, sizeof(T)*3*3);
01021 //
01022 //          //set a column of the matrix to i_b
01023 //          modifiedMatrix[0][column] = i_b[0];
01024 //          modifiedMatrix[1][column] = i_b[1];
01025 //          modifiedMatrix[2][column] = i_b[2];
01026 //
01027 //          o_x[column] = computeDeterminant(modifiedMatrix) / origDet;
01028 //        }
01029 //      }
01030 //      else {
01031 //        std::cerr << "Warning: Linear dependent eigenvectors! (using Jacobi Solver)" << std::endl;
01032 //        std::cerr << "Determinant: " << origDet << std::endl;
01033 //        std::cerr << i_matrix[0][0] << "\t" << i_matrix[0][1] << "\t" << i_matrix[0][2] << std::endl;
01034 //        std::cerr << i_matrix[1][0] << "\t" << i_matrix[1][1] << "\t" << i_matrix[1][2] << std::endl;
01035 //        std::cerr << i_matrix[2][0] << "\t" << i_matrix[2][1] << "\t" << i_matrix[2][2] << std::endl;
01036 //#if 0
01037 //        T xTemp[3];
01038 //        xTemp[0] = xTemp[1] = xTemp[2] = 0.;
01039 //        for (int m = 0; m < 20; m++) {
01040 //          for (int row = 0; row < 3; row++) {
01041 //            o_x[row] = 0.;
01042 //            for (int col = 0; col < 3; col++) {
01043 //              if (col != row)
01044 //                o_x[row] += i_matrix[row][col]*xTemp[col];
01045 //            }
01046 //            o_x[row] = (i_b[row] - o_x[row])/i_matrix[row][row];
01047 //          }
01048 //
01049 //          if (fabs(o_x[0]-xTemp[0]) + fabs(o_x[1]-xTemp[1]) + fabs(o_x[2]-xTemp[2]) < zeroTol*10.)
01050 //            break;
01051 //          else {
01052 //            std::cout << "Error: " << fabs(o_x[0]-xTemp[0]) + fabs(o_x[1]-xTemp[1]) + fabs(o_x[2]-xTemp[2]) << std::endl;
01053 //            xTemp[0] = o_x[0];
01054 //            xTemp[1] = o_x[1];
01055 //            xTemp[2] = o_x[2];
01056 //          }
01057 //        }
01058 //        std::cout << "Solution:" << std::endl;
01059 //        std::cout << "\t" << o_x[0] << std::endl;
01060 //        std::cout << "x=\t" << o_x[1] << std::endl;
01061 //        std::cout << "\t" << o_x[2] << std::endl;
01062 //        std::cout << "*********" << std::endl;
01063 //        std::cout << "\t" << i_b[0] << std::endl;
01064 //        std::cout << "b=\t" << i_b[1] << std::endl;
01065 //        std::cout << "\t" << i_b[2] << std::endl;
01066 //        std::cout << "***A*x****" << std::endl;
01067 //        for (int row = 0; row < 3; row++) {
01068 //          std::cout << "\t" <<
01069 //              i_matrix[row][0] * o_x[0] +
01070 //              i_matrix[row][1] * o_x[1] +
01071 //              i_matrix[row][2] * o_x[2] << std::endl;
01072 //        }
01073 //#endif
01074 //        assert(false);
01075 //      }
01076 //#endif
01077     }
01078 //#if 0
01079 //    /**
01080 //     * Compute the determinant of a 3x3 matrix
01081 //     * using formula: |A|=a11 * (a22a33-a32a23) + a12 * (a23a31-a33a21) + a13 * (a21a32-a31a22)
01082 //     *                |A|=a00 * (a11a22-a21a12) + a01 * (a12a20-a22a10) + a02 * (a10a21-a20a11)
01083 //     */
01084 //    T computeDeterminant(T matrix[3][3]) const
01085 //    {
01086 //
01087 //      T determinant =     matrix[0][0] * ( matrix[1][1] * matrix[2][2] - matrix[2][1] * matrix[1][2] )
01088 //                                          + matrix[0][1] * ( matrix[1][2] * matrix[2][0] - matrix[2][2] * matrix[1][0] )
01089 //                                          + matrix[0][2] * ( matrix[1][0] * matrix[2][1] - matrix[2][0] * matrix[1][1] );
01090 //
01091 //      return determinant;
01092 //    }
01093 //#endif
01094 
01095 
01096 #undef sqrt_g_hLeft
01097 #undef sqrt_g_hRight
01098 #undef sqrt_hLeft
01099 #undef sqrt_hRight
01100 
01101 };
01102 
01103 #undef complexsteadystatewave
01104 #undef correctrarefactions
01105 
01106 
01107 #endif /* AUGRIE_HPP_ */
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends