Chombo + EB: IndexedMomentsImplem.H Source File

00001 #ifdef CH_LANG_CC
00002 /*
00003  *      _______              __
00004  *     / ___/ /  ___  __ _  / /  ___
00005  *    / /__/ _ \/ _ \/  V \/ _ \/ _ \
00006  *    \___/_//_/\___/_/_/_/_.__/\___/
00007  *    Please refer to Copyright.txt, in Chombo's root directory.
00008  */
00009 #endif
00010 
00011 #ifndef _INDEXEDMOMENTSIMPLEM_H_
00012 #define _INDEXEDMOMENTSIMPLEM_H_
00013 
00014 #include "MomentIterator.H"
00015 #include "EB_TYPEDEFS.H"
00016 #include "IndexTM.H"
00017 #include "Factorial.H"
00018 #include "CH_Timer.H"
00019 #include "parstream.H"
00020 #include <map>
00021 #include <utility>
00022 #include "NamespaceHeader.H"
00023 
00024 template <int Dim, int P> bool                      IndexedMoments<Dim, P>::s_staticsSet     = false;
00025 template <int Dim, int P> Vector<IndexTM<int,Dim> > IndexedMoments<Dim, P>::s_multiIndicies  = Vector<IndexTM<int, Dim> >();
00026 template <int Dim, int P> int                       IndexedMoments<Dim, P>::s_size           = 0;
00027 using std::endl;
00028 
00029 ///
00030 /**
00031    Moments are centered at the center of the cell.
00032    For each of these moments I shift them to the lowest 
00033    corner of the cell, where I know what the bounds of the integrals 
00034    is (lower bound always zero, upper bound = dx^d dx^px dx^py dx^pz
00035    If the shifted moment is out of bounds, I bound it throw a message.
00036    tolerance is about verbosity.  
00037 **/
00038 template<int Dim, int ORDER>
00039 void
00040 checkMoments(IndexedMoments<Dim, ORDER> & a_moments,
00041              Vector<IndexTM<int, Dim> > & a_bogusPowers,
00042              const Real                 & a_dx,
00043              const Real                 & a_tolerance,
00044              const bool                 & a_ebMoment,
00045              const bool                 & a_bindMoments)
00046 {
00047   MomentIterator<Dim, ORDER> momit;
00048   //eb moments can mysteriously get the wrong signs
00049   if(a_ebMoment && a_bindMoments)
00050   {
00051     if(a_moments[IndexTM<int,Dim>::Zero] < 0.)
00052     {
00053       a_moments *= -1.0;
00054     }
00055   }
00056   a_bogusPowers.resize(0);
00057   Real cellvol = 1;
00058   for(int idir = 0; idir < Dim; idir++)
00059     {
00060       cellvol *= a_dx;
00061     }
00062 
00063   //shift moments to lowest corner 
00064   IndexTM<Real, Dim> shiftvec;
00065   shiftvec.setAll(0.5*a_dx);
00066   a_moments.shift(shiftvec);
00067 
00068   //because they are all shifted, to the lowest 
00069   //min is always zero.
00070   const Real zerval = 0;
00071   for(momit.reset(); momit.ok(); ++momit)
00072     {
00073       const IndexTM<int, Dim>& p = momit();
00074       //this makes max = dx^p
00075       Real maxval = POW<Dim>(a_dx, p);
00076       //because this is a vol integral...
00077       maxval *= cellvol;
00078       //this is a reference so it will change the input
00079       Real& momval = a_moments[p];
00080       if(momval < zerval)
00081         {
00082           if((zerval - momval) > a_tolerance*maxval)
00083             {
00084               a_bogusPowers.push_back(p);
00085               //              MayDay::Error("minval break low");
00086             }
00087           if(a_bindMoments)
00088             {
00089               momval = zerval;
00090             }
00091         }
00092       else if((momval > maxval) && (!a_ebMoment))
00093         {
00094           if((momval - maxval) > a_tolerance*maxval)
00095             {
00096               a_bogusPowers.push_back(p);
00097               //              MayDay::Error("maxval break");
00098             }
00099           if(a_bindMoments)
00100             {
00101               momval = maxval;
00102             }
00103         }
00104     }
00105 
00106   //shift moments back to highest corner 
00107   shiftvec.setAll(-a_dx);
00108   a_moments.shift(shiftvec);
00109 
00110   //because the moments are shifted to the highest corner
00111   //even powers will always be positive, odd ones will always be negtive
00112   for(momit.reset(); momit.ok(); ++momit)
00113     {
00114       const IndexTM<int, Dim>& p = momit();
00115       int psum = p.sum();
00116       bool pEven = (psum % 2 == 0);
00117       Real& momval = a_moments[p];
00118       //this makes max = dx^p
00119       Real maxval = POW<Dim>(a_dx, p);
00120       //because this is a vol integral...
00121       maxval *= cellvol;
00122       if(pEven && (momval < zerval) && (!a_ebMoment))
00123         {
00124           //p has an  sum and therefore the moment should be positive
00125           if(momval < -a_tolerance*maxval )
00126             {
00127               a_bogusPowers.push_back(p);
00128               //              MayDay::Error("high corner break even");
00129             }
00130           if(a_bindMoments)
00131             {
00132               momval = 0;
00133             }
00134         }
00135       else if(!pEven && (momval > zerval) && (!a_ebMoment))
00136         {
00137           //p has an odd sum and therefore the moment should be negative
00138           if(momval > a_tolerance*maxval)
00139             {
00140               a_bogusPowers.push_back(p);
00141               //              MayDay::Error("high corner break odd");
00142             }
00143           if(a_bindMoments)
00144             {
00145               momval = 0;
00146             }
00147         }
00148     }
00149 
00150   //shift moments back to cell center
00151   shiftvec.setAll(0.5*a_dx);
00152   a_moments.shift(shiftvec);
00153 }
00154 
00155 template <int Dim, int P>
00156 const int IndexedMoments<Dim, P>::s_max_sizes[/* D */][CH_IM_MAX_POWER+1] = 
00157   {
00158     {1, 2,  3,  4,  5,  6,  7,   8,   9}, // D=1
00159     {1, 3,  6, 10, 15, 21, 28,  36,  45}, // D=2
00160     {1, 4, 10, 20, 35, 56, 84, 120, 165}  // D=3
00161   }; 
00162 using std::endl; 
00163 /////
00164 template <const int Dim, const int P>
00165 Real
00166 IndexedMoments<Dim, P>:: 
00167 getMoment(const     IndexTM<int, Dim>        & a_mono,
00168           const map<IndexTM<int, Dim>,  Real>& a_mapin) const
00169 {
00170   /**/
00171   typename std::map<IndexTM<int, Dim>,  Real>::const_iterator iterinator;
00172   iterinator = a_mapin.find(a_mono);
00173 
00174   Real moment = 1.23456789e10;
00175 
00176   if (iterinator != a_mapin.end())
00177     {
00178       moment = iterinator->second;
00179 
00180     }
00181   else
00182     {
00183       pout() << "monomial = " ;
00184       for(int idir = 0; idir < Dim; idir++)
00185         {
00186           pout() << a_mono[idir] << " ";
00187         }
00188       pout() << endl;
00189       MayDay::Abort("getMoments in IndexedMoments: moment not found in map");
00190     }
00191   return moment;
00192 }
00193 
00194 ///multiply each component by factor
00195 template <int Dim, int P>
00196 IndexedMoments<Dim, P>&
00197 IndexedMoments<Dim, P>:: 
00198 operator*=(const Real& a_factor)
00199 {
00200   for(MomentIterator<Dim,P> it; it.ok(); ++it)
00201     {
00202       (*this)[it()] *= a_factor;
00203     }
00204   return *this;
00205 }
00206 ///add compoenentwise
00207 template <int Dim, int P>
00208 IndexedMoments<Dim, P>&
00209 IndexedMoments<Dim, P>:: 
00210 operator+=(const IndexedMoments<Dim, P>& a_input)
00211 {
00212   for(MomentIterator<Dim,P> it; it.ok(); ++it)
00213     {
00214       (*this)[it()] += a_input[it()];
00215     }
00216   return *this;
00217 }
00218 ///shift moment by the input distance.
00219 template <int Dim, int P>
00220 void
00221 IndexedMoments<Dim, P>:: 
00222 shift(const IndexTM<Real, Dim>& a_distance)
00223 {
00224   CH_TIME("IndexedMoment::shift");
00225   IndexTM<Real, Dim> x = a_distance;
00226 
00227   IndexedMoments<Dim, P> retval = *this;
00228   for(MomentIterator<Dim,P> itouter; itouter.ok(); ++itouter)
00229     {
00230       IndexTM<int, Dim> r = itouter();
00231       // - use the binomial shift rule to move them to a_x0
00232       Real moment = 0.0;
00233       for(MomentIterator<Dim,P> itinner; itinner.ok(); ++itinner)
00234         {
00235           IndexTM<int, Dim> q = itinner();
00236          
00237           if(q.componentwiseLE(r))
00238             {
00239               IndexTM<int, Dim> p = r - q;
00240               Real m0 = (*this)[q];
00241               Real xpow = 1;
00242               for(int idir = 0; idir < Dim; idir++)
00243                 {
00244                   xpow *=  POW(x[idir], p[idir]);
00245                 }
00246               
00247               moment += pCk<Dim>(r, q) * xpow * m0;
00248             
00249             } //end if (q <= r)
00250         }//end loop over inner moments
00251       retval[r] = moment;
00252     }
00253   *this = retval;
00254 }
00255 
00256 /*******/
00257 /// for use with irregnode
00258 template <int Dim, int P>
00259 IndexedMoments<Dim, P>& 
00260 IndexedMoments<Dim, P>:: 
00261 operator=(const map<IndexTM<int, Dim>,  Real>& a_mapin)
00262 {
00263   
00264   for(MomentIterator<Dim, P> iter; iter.ok(); ++iter)
00265     {
00266       IndexTM<int, Dim> mono = iter();
00267       Real moment= getMoment(iter(), a_mapin);
00268       (*this)[iter()] = moment;
00269     }
00270   return *this;
00271 }
00272 
00273 template <int Dim, int P>
00274 void
00275 IndexedMoments<Dim, P>:: 
00276 spout() const
00277 {
00278   pout() << "index \t  moment  " << endl;
00279   for(MomentIterator<Dim, P> iter; iter.ok(); ++iter)
00280     {
00281       pout() <<  iter() << "\t" << (*this)[iter()] << endl;
00282     }
00283 }
00284 
00285 template <int Dim>
00286 bool allPositive(const IndexTM<int, Dim>& a_index) 
00287 {
00288   bool allpos = true;
00289   for(int idir = 0; idir < Dim; idir++)
00290     {
00291       if(a_index[idir] < 0)
00292         {
00293           allpos = false;
00294         }
00295     }
00296   return allpos;
00297 }
00298 
00299 
00300 template <int Dim, int P>
00301 int
00302 IndexedMoments<Dim, P>::
00303 indexOf(IndexTM<int,Dim> a_index)
00304 {
00305   CH_assert(a_index.sum() <= P);
00306   CH_assert(allPositive(a_index));
00307 
00308   // Formula for the layout index with 0th dim first, then 1st, etc.
00309   int index = a_index[0];
00310   int rem = P;
00311   for (int d=Dim-1; d > 0; --d)
00312     {
00313       index += s_max_sizes[d][rem] - s_max_sizes[d][rem-a_index[d]];
00314       rem -= a_index[d];
00315     }
00316 
00317   return index;
00318 }
00319 /*******/
00320 template <int Dim, int P>
00321 void
00322 IndexedMoments<Dim, P>::
00323 setStatics() 
00324 {
00325   setSize();
00326   setMultiIndicies();
00327   s_staticsSet= true;
00328 }
00329 
00330 /*******/
00331 template <int Dim, int P>
00332 IndexedMoments<Dim, P>::
00333 IndexedMoments() 
00334 {
00335   if(!s_staticsSet)
00336     {
00337       setStatics();
00338     }
00339 
00340   m_moms.resize(s_size);
00341   m_isRegular = 0;
00342 }
00343 
00344 /*******/
00345 template <int Dim, int P>
00346 void
00347 IndexedMoments<Dim, P>::
00348 setRegular(const Real a_dx)
00349 {
00350   CH_assert(a_dx > 0);
00351   m_isRegular = 1;
00352 
00353   // Get the set of multi-indexes in the correct order
00354   const Vector<IndexTM<int,Dim> >& mix = IndexedMoments<Dim, P>::getMonomialPowers();
00355 
00356   // pout() << "Regular moments:" << endl;
00357   for (int ix=0; ix < mix.size(); ++ix)
00358     {
00359       IndexTM<int,Dim> index = mix[ix];
00360 
00361       bool even = true;
00362       for (int d=0; (even) && (d < Dim); ++d)
00363         even = even && !(index[d] % 2);
00364 
00365       //! Only moments with an even multi-index are non-zero
00366       Real moment = 0;
00367       if (even)
00368         {
00369           moment = POW(0.5 * a_dx, index.sum()) * POW(a_dx, Dim);
00370           for (int d=0; d < Dim; ++d)
00371             moment /= (Real) (index[d]+1);
00372         }
00373       m_moms[ix] = moment;
00374 
00375       // pout() << ix << ": " << index << " = " << moment << endl;
00376     }
00377 
00378 }
00379 /*******/
00380 template <int Dim, int P>
00381 void
00382 IndexedMoments<Dim, P>::
00383 setSize()
00384 {
00385   CH_assert(Dim <= SpaceDim);
00386   CH_assert(Dim > 0);
00387   CH_assert(P <  CH_IM_MAX_POWER);
00388   CH_assert(P >= 0);
00389   
00390 
00391   s_size = s_max_sizes[Dim-1][P];
00392 }
00393 /*******/
00394 template <int Dim, int P>
00395 void
00396 IndexedMoments<Dim, P>::
00397 setMultiIndicies()
00398 {
00399   s_multiIndicies.resize(s_size);
00400 
00401   IndexTM<int,Dim> index = IndexTM<int,Dim>::Zero;
00402   for (int ix=0; ix < s_size; ++ix)
00403     {
00404       // If the sum is too large, shift extras to the right
00405       for (int d=0; (index.sum() > P) && (d < Dim-1); ++d)
00406         {
00407           index[d] = 0;
00408           ++index[d+1];
00409         }
00410       s_multiIndicies[ix] = index;
00411       ++index[0];
00412     }
00413 }
00414 /****/
00415 template <int Dim, int P>
00416 void 
00417 IndexedMoments<Dim, P>::
00418 setToTruncatedMultiply(const IndexedMoments<Dim, P> & a_CA,
00419                        const IndexedMoments<Dim, P> & a_CB)
00420 {
00421   this->setToZero();
00422   //here we are computing sum(m_p) = integral(ni(x-xbar)^p1/p! (x-xbar)^p2)
00423   // = sum_p2 ((n^p2)/p2!) m_p+p12
00424   for(MomentIterator<Dim, P> momitOuter; momitOuter.ok(); ++momitOuter)
00425     {
00426       //this is where the answer goes 
00427       IndexTM<int,SpaceDim>  po = momitOuter();
00428 
00429       for(MomentIterator<Dim, P> momitInner; momitInner.ok(); ++momitInner)
00430         {
00431           IndexTM<int,SpaceDim> pi = momitInner();
00432           IndexTM<int,SpaceDim> psum = po + pi;
00433           if(psum.sum() <= P)
00434             {
00435               Real incr = a_CA[pi]*a_CB[po];
00436               (*this)[psum] += incr;
00437             }
00438         }
00439     }
00440 
00441 }
00442 /****/
00443 template <int Dim, int P>
00444 void 
00445 IndexedMoments<Dim, P>::
00446 divideByFactorial()
00447 {
00448   for(MomentIterator<Dim, P> momit; momit.ok(); ++momit)
00449     {
00450       IndexTM<int,SpaceDim>  p = momit();
00451       Real factorial = pfactorial(p);
00452       (*this)[p] /= factorial;
00453     }
00454 }
00455 /****/
00456 template <int Dim, int P>
00457 void 
00458 IndexedMoments<Dim, P>::
00459 multiplyByFactorial()
00460 {
00461   for(MomentIterator<Dim, P> momit; momit.ok(); ++momit)
00462     {
00463       IndexTM<int,SpaceDim>  p = momit();
00464       Real factorial = pfactorial(p);
00465       (*this)[p] *= factorial;
00466     }
00467 }
00468 /*******/
00469 #include "NamespaceFooter.H"
00470 #endif