Chombo + EB: IndexTMI.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 _INDEXTMI_H_
00012 #define _INDEXTMI_H_
00013 
00014 #include "SPACE.H"
00015 #include <iostream>
00016 using std::ostream;
00017 using std::istream;
00018 using std::ws;
00019 
00020 #include "MayDay.H"
00021 #include "Misc.H"
00022 #include "IndexTM.H"
00023 #include "parstream.H"
00024 #include <cmath>
00025 #include <algorithm>
00026 #include <functional>
00027 #include <numeric>
00028 #include "Metaprograms.H"
00029 #include "NamespaceHeader.H"
00030 
00031 
00032 template<typename T, int N> ostream&
00033 operator<< (ostream&       os,
00034             const IndexTM<T,N>& p)
00035 {
00036     os << '(' << p[0];
00037     for( int i=1;i<N;++i )
00038     {
00039         os << ',' << p[i];
00040     }
00041     os << ')';
00042     if (os.fail())
00043         MayDay::Error("operator<<(ostream&,Index&) failed");
00044     return os;
00045 }
00046 
00047 //
00048 // Copied from <Utility.H>
00049 //
00050 #define CH_IGNORE_MAX 100000
00051 
00052 template<typename T, int N> istream&
00053 operator>> (istream& is,
00054             IndexTM<T,N>& p)
00055 {
00056     is >> ws;
00057     char c;
00058     is >> c;
00059     is.putback(c);
00060     if (c == '(')
00061     {
00062         is.ignore(CH_IGNORE_MAX, '(') >> p[0];
00063         for( int i=1;i<N;++i )
00064         {
00065             is.ignore(CH_IGNORE_MAX, ',') >> p[i];
00066         }
00067         is.ignore(CH_IGNORE_MAX, ')');
00068     }
00069     else if (c == '<')
00070     {
00071         is.ignore(CH_IGNORE_MAX, '<') >> p[0];
00072         for( int i=1;i<N;++i )
00073         {
00074             is.ignore(CH_IGNORE_MAX, ',') >> p[1];
00075         }
00076         is.ignore(CH_IGNORE_MAX, '>');
00077     }
00078     else
00079         MayDay::Error("operator>>(istream&,Index&): expected \'(\' or \'<\'");
00080 
00081     if (is.fail())
00082         MayDay::Error("operator>>(istream&,Index&) failed");
00083 
00084     return is;
00085 }
00086 
00087 
00088 template<typename T, int N> void
00089 IndexTM<T,N>::printOn (ostream& os) const
00090 {
00091     os << "Index: " << *this << '\n';
00092 }
00093 
00094 template<typename T, int N> void
00095 IndexTM<T,N>::p() const
00096 {
00097     pout() << *this << '\n';
00098 }
00099 
00100 template<typename T, int N> void
00101 IndexTM<T,N>::dumpOn (ostream& os) const
00102 {
00103     os << "Index " << *this << '\n';
00104 }
00105 
00106 //
00107 // Static object initialization.
00108 //
00109 template<typename T, int N> int
00110 IndexTM<T,N>::InitStatics()
00111 {
00112   IndexTM<T,N>& pz = const_cast<IndexTM<T,N>&>( IndexTM<T,N>::Zero );
00113 
00114   pz.setAll(0);
00115 
00116   IndexTM<T,N>& pu = const_cast<IndexTM<T,N>&>( IndexTM<T,N>::Unit );
00117 
00118   pu.setAll(1);
00119 
00120   // No danger of Index::Zero and Unit not having been allocated, as ARM section
00121   // 3.4 says "The initialization of nonlocal static objects in a translation unit
00122   // is done before the first use of any function...defined in that translation
00123   // unit."
00124   //
00125   // Had to go through the const_cast stuff because it's nice to be able to declare
00126   // Index::Zero and Index::Unit as const.
00127 
00128   return 0; // arbitrary
00129 }
00130 
00131 //
00132 // Inlines.
00133 //
00134 template<typename T, int N> inline
00135 IndexTM<T,N>::IndexTM(const T *a)
00136   : GenericArithmeticable< T, IndexTM<T,N> >( this )
00137 {
00138   memcpy( vect, a, N*sizeof(T) );
00139 }
00140 
00141 template<typename T, int N> inline
00142 IndexTM<T,N>::IndexTM(const IndexTM<T,N> &iv)
00143   : GenericArithmeticable< T, IndexTM<T,N> >( this )
00144 {
00145   memcpy(vect, iv.vect, N*sizeof(T) );
00146 }
00147 
00148 template<typename T, int N> inline
00149 IndexTM<T,N>&
00150 IndexTM<T,N>::operator=(const IndexTM<T,N> &iv)
00151 {
00152   memcpy( vect, iv.vect, N*sizeof(T) );
00153   return *this;
00154 }
00155 
00156 template<typename T, int N> inline
00157 T&
00158 IndexTM<T,N>::operator[] (int i)
00159 {
00160   CH_assert(i>=0 && i < N);
00161   return vect[i];
00162 }
00163 
00164 template<typename T, int N> inline
00165 T
00166 IndexTM<T,N>::operator[] (int i) const
00167 {
00168   CH_assert(i>=0 && i < N);
00169   return vect[i];
00170 }
00171 
00172 template<typename T, int N> inline
00173 void
00174 IndexTM<T,N>::setVal (int i,
00175                       T val)
00176 {
00177   CH_assert(i >=0 && i < N);
00178   vect[i] = val;
00179 }
00180 
00181 template<typename T, int N> inline
00182 void
00183 IndexTM<T,N>::setAll (T val)
00184 {
00185   for(int i=0; i<N; i++)
00186     vect[i] = val;
00187 }
00188 
00189 
00190 template<typename T, int N> inline
00191 const T*
00192 IndexTM<T,N>::dataPtr() const
00193 {
00194   return vect;
00195 }
00196 
00197 template<typename T, int N> inline
00198 T*
00199 IndexTM<T,N>::dataPtr()
00200 {
00201   return vect;
00202 }
00203 
00204 template<typename T, int N> inline
00205 const T*
00206 IndexTM<T,N>::getVect () const
00207 {
00208   return vect;
00209 }
00210 
00211 template<typename T, int N> inline
00212 bool
00213 IndexTM<T,N>::operator== (const IndexTM& p) const
00214 {
00215   return Metaprograms::pointwiseCompare<N,T,std::equal_to<T> >(vect,p.vect);
00216 }
00217 
00218 template<typename T, int N> inline
00219 bool
00220 IndexTM<T,N>::operator!= (const IndexTM& p) const
00221 {
00222   return ! operator==(p);
00223 }
00224 
00225 template<typename T, int N> inline
00226 bool
00227 IndexTM<T,N>::lexLT (const IndexTM &s) const
00228 {
00229   return Metaprograms::LexLT<N,T>()(vect,s.vect);
00230 }
00231 
00232 template<typename T, int N> inline
00233 bool
00234 IndexTM<T,N>::lexGT (const IndexTM& s) const
00235 {
00236   return ! Metaprograms::LexLT<N,T>()(vect,s.vect);
00237 }
00238 
00239 template<typename T, int N> inline
00240 IndexTM<T,N>
00241 IndexTM<T,N>::operator+ () const
00242 {
00243   return IndexTM<T,N>(*this);
00244 }
00245 
00246 template<typename T, int N> inline
00247 IndexTM<T,N>
00248 IndexTM<T,N>::operator- () const
00249 {
00250   IndexTM<T,N> result(*this);
00251   for( int i=0;i<N;++i )
00252     {
00253       result.vect[i] *= -1;
00254     }
00255   return result;
00256 }
00257 
00258 template<typename T, int N> inline
00259 T
00260 IndexTM<T,N>::dotProduct ( const IndexTM<T,N>& a_rhs ) const
00261 {
00262   return Metaprograms::InnerProduct<N,T,T,
00263                                     std::plus<T>,
00264                                     std::multiplies<T> >()( vect,
00265                                                             a_rhs.vect );
00266 }
00267 
00268 template<typename T, int N> inline
00269 T
00270 IndexTM<T,N>::sum () const
00271 {
00272   return Metaprograms::Accum<N,T,std::plus<T> >()( vect );
00273 }
00274 
00275 template<typename T, int N> inline
00276 T
00277 IndexTM<T,N>::product () const
00278 {
00279   return Metaprograms::Accum<N,T,std::multiplies<T> >()( vect );
00280 }
00281 
00282 template<typename T, int N> template<typename OP> bool
00283 IndexTM<T,N>::operatorCompare(const IndexTM<T,N>& p, const OP& op) const
00284 {
00285   return Metaprograms::pointwiseCompare<N,T,OP>(vect,p.vect);
00286 }
00287 
00288 
00289 template<typename T, int N> template<typename OP> inline
00290 IndexTM<T,N>&
00291 IndexTM<T,N>::operatorOpEquals(const T& s, const OP& op)
00292 {
00293   Metaprograms::Transform<N,T,OP>()(vect,s);
00294   return *this;
00295 }
00296 
00297 template<typename T, int N> template<typename OP> inline
00298 IndexTM<T,N>&
00299 IndexTM<T,N>::operatorOpEquals(const IndexTM<T,N>& p, const OP& op)
00300 {
00301   Metaprograms::Transform<N,T,OP>()( vect, p.vect );
00302   return *this;
00303 }
00304 
00305 template<typename T, int N> inline
00306 IndexTM<T,N>&
00307 IndexTM<T,N>::reciprocal()
00308 {
00309   std::transform( vect, vect+N, vect, bind1st(std::divides<T>(),T(1)) );
00310   return *this;
00311 }
00312 
00313 
00314 template<typename T> static bool abscompare(const T& a, const T& b)
00315 {
00316   return ::fabs(a) < ::fabs(b);
00317 }
00318 template<typename T, int N> inline
00319 int
00320 IndexTM<T,N>::minDir (bool a_doAbs ) const
00321 {
00322   if( a_doAbs )
00323     {
00324       return std::min_element( vect, vect+N, std::ptr_fun(abscompare<T>) ) - vect;
00325     }
00326   else
00327     {
00328       return std::min_element( vect, vect+N ) - vect;
00329     }
00330 }
00331 template<typename T, int N> inline
00332 int
00333 IndexTM<T,N>::maxDir (bool a_doAbs ) const
00334 {
00335   if( a_doAbs )
00336     {
00337       return std::max_element( vect, vect+N, std::ptr_fun(abscompare<T>) ) - vect;
00338     }
00339   else
00340     {
00341       return std::max_element( vect, vect+N ) - vect;
00342     }
00343 }
00344 
00345 template<typename T> static T ourmin(T a, T b) { return a<b ? a:b; }
00346 template<typename T> static T ourmax(T a, T b) { return a>b ? a:b; }
00347 
00348 template<typename T, int N> inline
00349 IndexTM<T,N>&
00350 IndexTM<T,N>::min (const IndexTM<T,N>& p)
00351 {
00352   std::transform( vect, vect+N, p.vect, vect,
00353                   std::ptr_fun(ourmin<T>) );
00354   return *this;
00355 }
00356 
00357 template<typename T, int N> inline
00358 IndexTM<T,N>&
00359 IndexTM<T,N>::max (const IndexTM<T,N>& p)
00360 {
00361   std::transform( vect, vect+N, p.vect, vect,
00362                   std::ptr_fun(ourmax<T>) );
00363   return *this;
00364 }
00365 
00366 
00367 template<typename T, int N> inline
00368 IndexTM<T,N>&
00369 IndexTM<T,N>::scale (T s)
00370 {
00371     return (*this) *= s;
00372 }
00373 
00374 template<typename T, int N> inline
00375 IndexTM<T,N>&
00376 IndexTM<T,N>::reflect (T ref_ix,
00377                        int idir)
00378 {
00379   CH_assert(idir >= 0 && idir < N);
00380   vect[idir] = -vect[idir] + 2*ref_ix;
00381   return *this;
00382 }
00383 
00384 template<typename T, int N> inline
00385 IndexTM<T,N>&
00386 IndexTM<T,N>::shift (int coord,
00387                   T s)
00388 {
00389   CH_assert(coord >= 0 && coord < N);
00390   vect[coord] += s;
00391   return *this;
00392 }
00393 
00394 template<typename T, int N> inline
00395 IndexTM<T,N>&
00396 IndexTM<T,N>::shift (const IndexTM<T,N>& iv)
00397 {
00398   return (*this) += iv;
00399 }
00400 
00401 template<typename T, int N> inline
00402 IndexTM<T,N>&
00403 IndexTM<T,N>::diagShift (T s)
00404 {
00405   return (*this) += s;
00406 }
00407 
00408 
00409 template<typename T, int N> inline
00410 IndexTM<T,N>
00411 scale (const IndexTM<T,N>& p,
00412        T                   s)
00413 {
00414   return p * s;
00415 }
00416 
00417 template<typename T, int N> inline
00418 IndexTM<T,N>
00419 diagShift(const IndexTM<T,N> &p, T s)
00420 {
00421   return p + s;
00422 }
00423 
00424 template<typename T, int N> inline
00425 IndexTM<T,N>
00426 min(const IndexTM<T,N>& p1,
00427     const IndexTM<T,N>& p2)
00428 {
00429   IndexTM<T,N> result( p1 );
00430   return result.min(p2);
00431 }
00432 
00433 template<typename T, int N> inline
00434 IndexTM<T,N>
00435 max(const IndexTM<T,N>& p1,
00436     const IndexTM<T,N>& p2)
00437 {
00438   IndexTM<T,N> result( p1 );
00439   return result.max(p2);
00440 }
00441 
00442 template<typename T, int N> inline
00443 IndexTM<T,N>
00444 BASISV_TM(int dir)
00445 {
00446   CH_assert(dir >= 0 && dir < N);
00447   IndexTM<T,N> tmp = IndexTM<T,N>::Zero ;
00448   tmp.dataPtr()[dir] = T(1);
00449   return tmp;
00450 }
00451 
00452 template<typename T, int N> inline
00453 IndexTM<T,N>
00454 reflect(const IndexTM<T,N>& a,
00455         T            ref_ix,
00456         int            idir)
00457 {
00458   IndexTM<T,N> result( a );
00459   return result.reflect( ref_ix, idir );
00460 }
00461 
00462 template<typename T> static T ourcoarsen( T a, T b )
00463 {
00464     return a<0 ? T(-::fabs(a+1))/b-1 : a/b;
00465 }
00466 
00467 template<typename T, int N> inline
00468 IndexTM<T,N>
00469 coarsen(const IndexTM<T,N>& p,
00470         T                   s)
00471 {
00472   IndexTM<T,N> result( p );
00473   return result.coarsen( s );
00474 }
00475 
00476 template<typename T, int N> inline
00477 IndexTM<T,N>
00478 coarsen(const IndexTM<T,N>& p1,
00479         const IndexTM<T,N>& p2)
00480 {
00481   IndexTM<T,N> result( p1 );
00482   return result.coarsen( p2 );
00483 }
00484 
00485 /*
00486 template<typename T, int N>
00487 IndexTM<T,N>::operator IndexTM<Real,N>()
00488 {
00489   IndexTM<Real,N> result;
00490   for( int i=0;i<N;++i ) result.dataPtr()[i] = Real(vect[i]);
00491   return result;
00492 }
00493 */
00494 
00495 template<typename T, int N> inline
00496 IndexTM<T,N>&
00497 IndexTM<T,N>::coarsen (T s)
00498 {
00499   CH_assert(s > 0);
00500   std::transform( vect, vect+N, vect,
00501                   std::bind2nd(std::ptr_fun(ourcoarsen<T>),s) );
00502   return *this;
00503 }
00504 
00505 template<typename T, int N> inline
00506 IndexTM<T,N>&
00507 IndexTM<T,N>::coarsen (const IndexTM<T,N>& p)
00508 {
00509   CH_assert(p > (IndexTM<T,N>::Zero));
00510   std::transform( vect, vect+N, p.vect, vect, std::ptr_fun(ourcoarsen<T>) );
00511   return *this;
00512 }
00513 
00514 template<typename T, int N>
00515 void
00516 IndexTM<T,N>::linearIn( const void* inBuf )
00517 {
00518   memcpy( vect, (T*)inBuf, N*sizeof(T) );
00519 }
00520 template<typename T, int N>
00521 void
00522 IndexTM<T,N>::linearOut( void* inBuf ) const
00523 {
00524   memcpy( (T*)inBuf, vect, N*sizeof(T) );
00525 }
00526 
00527 
00528 template<typename T, int N>
00529 IndexTM<T,N>::IndexTM( T i )
00530   : GenericArithmeticable< T, IndexTM<T,N> >( this )
00531 {
00532   STATIC_ASSERT( N==1 );
00533   vect[0] = i;
00534 }
00535 template<typename T, int N>
00536 IndexTM<T,N>::IndexTM( T i, T j )
00537   : GenericArithmeticable< T, IndexTM<T,N> >( this )
00538 {
00539   STATIC_ASSERT( N==2 );
00540   vect[0] = i;
00541   vect[1] = j;
00542 }
00543 template<typename T, int N>
00544 IndexTM<T,N>::IndexTM( T i, T j, T k )
00545   : GenericArithmeticable< T, IndexTM<T,N> >( this )
00546 {
00547   STATIC_ASSERT( N==3 );
00548   vect[0] = i;
00549   vect[1] = j;
00550   vect[2] = k;
00551 }
00552 template<typename T, int N>
00553 IndexTM<T,N>::IndexTM( T i, T j, T k, T m )
00554   : GenericArithmeticable< T, IndexTM<T,N> >( this )
00555 {
00556   STATIC_ASSERT( N==4 );
00557   vect[0] = i;
00558   vect[1] = j;
00559   vect[2] = k;
00560   vect[3] = m;
00561 }
00562 template<typename T, int N>
00563 IndexTM<T,N>::IndexTM( T i, T j, T k, T m, T n )
00564   : GenericArithmeticable< T, IndexTM<T,N> >( this )
00565 {
00566   STATIC_ASSERT( N==5 );
00567   vect[0] = i;
00568   vect[1] = j;
00569   vect[2] = k;
00570   vect[3] = m;
00571   vect[4] = n;
00572 }
00573 template<typename T, int N>
00574 IndexTM<T,N>::IndexTM( T i, T j, T k, T m, T n, T p )
00575   : GenericArithmeticable< T, IndexTM<T,N> >( this )
00576 {
00577   STATIC_ASSERT( N==6 );
00578   vect[0] = i;
00579   vect[1] = j;
00580   vect[2] = k;
00581   vect[3] = m;
00582   vect[4] = n;
00583   vect[5] = p;
00584 }
00585 
00586 
00587 template<typename T, int N> const IndexTM<T,N> IndexTM<T,N>::Zero;
00588 template<typename T, int N> const IndexTM<T,N> IndexTM<T,N>::Unit;
00589 static int s_dummyForIntVectCpp1( IndexTM<int,1>::InitStatics() );
00590 static int s_dummyForIntVectCpp2( IndexTM<int,2>::InitStatics() );
00591 static int s_dummyForIntVectCpp3( IndexTM<int,3>::InitStatics() );
00592 static int s_dummyForIntVectCpp4( IndexTM<int,4>::InitStatics() );
00593 static int s_dummyForIntVectCpp5( IndexTM<int,5>::InitStatics() );
00594 static int s_dummyForIntVectCpp6( IndexTM<int,6>::InitStatics() );
00595 static int s_dummyForRealVectCpp1( IndexTM<Real,1>::InitStatics() );
00596 static int s_dummyForRealVectCpp2( IndexTM<Real,2>::InitStatics() );
00597 static int s_dummyForRealVectCpp3( IndexTM<Real,3>::InitStatics() );
00598 static int s_dummyForRealVectCpp4( IndexTM<Real,4>::InitStatics() );
00599 static int s_dummyForRealVectCpp5( IndexTM<Real,5>::InitStatics() );
00600 static int s_dummyForRealVectCpp6( IndexTM<Real,6>::InitStatics() );
00601 // If Index::Zero and Index::Unit were pointers, we wouldn't need this extra
00602 // static int.  But they're objects, and the danger is that the initializations
00603 // right above here ("Index Index::Zero;" and "Index Index::Unit;") are hit
00604 // after the last call to Index::InitStatics, and in that case the
00605 // Index::Index() constructor could redefine Zero and Unit.  In fact, the way
00606 // things stand now, nothing bad would happen, because the Index::Index()
00607 // constructor doesn't assign anything to any of the data members.  But we don't
00608 // want to count on that always being the case.
00609 #include "NamespaceFooter.H"
00610 
00611 #endif // include guard