Intrepid
http://trilinos.sandia.gov/packages/docs/r11.2/packages/intrepid/src/Shared/MiniTensor/Intrepid_MiniTensor_Tensor3.t.h
00001 // @HEADER
00002 // ************************************************************************
00003 //
00004 //                           Intrepid Package
00005 //                 Copyright (2007) Sandia Corporation
00006 //
00007 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00008 // license for use of this work by or on behalf of the U.S. Government.
00009 //
00010 // Redistribution and use in source and binary forms, with or without
00011 // modification, are permitted provided that the following conditions are
00012 // met:
00013 //
00014 // 1. Redistributions of source code must retain the above copyright
00015 // notice, this list of conditions and the following disclaimer.
00016 //
00017 // 2. Redistributions in binary form must reproduce the above copyright
00018 // notice, this list of conditions and the following disclaimer in the
00019 // documentation and/or other materials provided with the distribution.
00020 //
00021 // 3. Neither the name of the Corporation nor the names of the
00022 // contributors may be used to endorse or promote products derived from
00023 // this software without specific prior written permission.
00024 //
00025 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00026 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00027 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00028 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00029 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00030 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00031 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00032 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00033 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00034 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00035 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00036 //
00037 // Questions: Alejandro Mota (amota@sandia.gov)
00038 //
00039 // ************************************************************************
00040 // @HEADER
00041 
00042 #if !defined(Intrepid_MiniTensor_Tensor3_t_h)
00043 #define Intrepid_MiniTensor_Tensor3_t_h
00044 
00045 namespace Intrepid {
00046 
00047   //
00048   // set dimension
00049   //
00050   //
00051   template<typename T>
00052   void
00053   Tensor3<T>::set_dimension(Index const N)
00054   {
00055     if (N == dimension) return;
00056 
00057     if (e != NULL) {
00058       delete [] e;
00059     }
00060 
00061     Index const
00062     number_components = N * N * N;
00063 
00064     e = new T[number_components];
00065 
00066     dimension = N;
00067 
00068     return;
00069   }
00070 
00071   //
00072   // 3rd-order tensor default constructor
00073   //
00074   template<typename T>
00075   Tensor3<T>::Tensor3() :
00076     dimension(0),
00077     e(NULL)
00078   {
00079     return;
00080   }
00081 
00082   //
00083   // 3rd-order tensor constructor with NaNs
00084   //
00085   template<typename T>
00086   Tensor3<T>::Tensor3(Index const N) :
00087     dimension(0),
00088     e(NULL)
00089   {
00090     set_dimension(N);
00091 
00092     Index const
00093     number_components = N * N * N;
00094 
00095     for (Index i = 0; i < number_components; ++i) {
00096       e[i] = not_a_number<T>();
00097     }
00098 
00099     return;
00100   }
00101 
00102   //
00103   // R^N 3rd-order tensor constructor with a scalar
00104   // \param s all components set to this scalar
00105   //
00106   template<typename T>
00107   Tensor3<T>::Tensor3(Index const N, T const & s) :
00108     dimension(0),
00109     e(NULL)
00110   {
00111     set_dimension(N);
00112 
00113     Index const
00114     number_components = N * N * N;
00115 
00116     for (Index i = 0; i < number_components; ++i) {
00117       e[i] = s;
00118     }
00119 
00120     return;
00121   }
00122 
00123   //
00124   // R^N copy constructor
00125   // 3rd-order tensor constructor from 3rd-order tensor
00126   // \param A from which components are copied
00127   //
00128   template<typename T>
00129   Tensor3<T>::Tensor3(Tensor3<T> const & A) :
00130     dimension(0),
00131     e(NULL)
00132   {
00133     Index const
00134     N = A.get_dimension();
00135 
00136     set_dimension(N);
00137 
00138     Index const
00139     number_components = N * N * N;
00140 
00141     for (Index i = 0; i < number_components; ++i) {
00142       e[i] = A.e[i];
00143     }
00144 
00145     return;
00146   }
00147 
00148   //
00149   // R^N 3rd-order tensor simple destructor
00150   //
00151   template<typename T>
00152   Tensor3<T>::~Tensor3()
00153   {
00154     if (e != NULL) {
00155       delete [] e;
00156     }
00157     return;
00158   }
00159 
00160   //
00161   // R^N 3rd-order tensor copy assignment
00162   //
00163   template<typename T>
00164   Tensor3<T> &
00165   Tensor3<T>::operator=(Tensor3<T> const & A)
00166   {
00167     if (this != &A) return *this;
00168 
00169     Index const
00170     N = A.get_dimension();
00171 
00172     set_dimension(N);
00173 
00174     Index const
00175     number_components = N * N * N;
00176 
00177     for (Index i = 0; i < number_components; ++i) {
00178       e[i] = A.e[i];
00179     }
00180 
00181     return *this;
00182   }
00183 
00184   //
00185   // 3rd-order tensor increment
00186   // \param A added to this tensor
00187   //
00188   template<typename T>
00189   Tensor3<T> &
00190   Tensor3<T>::operator+=(Tensor3<T> const & A)
00191   {
00192     Index const
00193     N = get_dimension();
00194 
00195     assert(A.get_dimension() == N);
00196 
00197     Index const
00198     number_components = N * N * N;
00199 
00200     for (Index i = 0; i < number_components; ++i) {
00201       e[i] += A.e[i];
00202     }
00203 
00204     return *this;
00205   }
00206 
00207   //
00208   // 3rd-order tensor decrement
00209   // \param A substracted from this tensor
00210   //
00211   template<typename T>
00212   Tensor3<T> &
00213   Tensor3<T>::operator-=(Tensor3<T> const & A)
00214   {
00215     Index const
00216     N = get_dimension();
00217 
00218     assert(A.get_dimension() == N);
00219 
00220     Index const
00221     number_components = N * N * N;
00222 
00223     for (Index i = 0; i < number_components; ++i) {
00224       e[i] -= A.e[i];
00225     }
00226 
00227     return *this;
00228   }
00229 
00230   //
00231   // R^N fill 3rd-order tensor with zeros
00232   //
00233   template<typename T>
00234   void
00235   Tensor3<T>::clear()
00236   {
00237     Index const
00238     N = get_dimension();
00239 
00240     Index const
00241     number_components = N * N * N;
00242 
00243     for (Index i = 0; i < number_components; ++i) {
00244       e[i] = 0.0;;
00245     }
00246 
00247     return;
00248   }
00249 
00250   //
00251   // 3rd-order tensor addition
00252   // \param A 3rd-order tensor
00253   // \param B 3rd-order tensor
00254   // \return \f$ A + B \f$
00255   //
00256   template<typename T>
00257   Tensor3<T>
00258   operator+(Tensor3<T> const & A, Tensor3<T> const & B)
00259   {
00260     Index const
00261     N = A.get_dimension();
00262 
00263     assert(B.get_dimension() == N);
00264 
00265     Tensor3<T>
00266     S(N);
00267 
00268     for (Index i = 0; i < N; ++i) {
00269       for (Index j = 0; j < N; ++j) {
00270         for (Index k = 0; k < N; ++k) {
00271           S(i,j,k) = A(i,j,k) + B(i,j,k);
00272         }
00273       }
00274     }
00275 
00276     return S;
00277   }
00278 
00279   //
00280   // 3rd-order tensor substraction
00281   // \param A 3rd-order tensor
00282   // \param B 3rd-order tensor
00283   // \return \f$ A - B \f$
00284   //
00285   template<typename T>
00286   Tensor3<T>
00287   operator-(Tensor3<T> const & A, Tensor3<T> const & B)
00288   {
00289     Index const
00290     N = A.get_dimension();
00291 
00292     assert(B.get_dimension() == N);
00293 
00294     Tensor3<T>
00295     S(N);
00296 
00297     for (Index i = 0; i < N; ++i) {
00298       for (Index j = 0; j < N; ++j) {
00299         for (Index k = 0; k < N; ++k) {
00300           S(i,j,k) = A(i,j,k) - B(i,j,k);
00301         }
00302       }
00303     }
00304 
00305     return S;
00306   }
00307 
00308   //
00309   // 3rd-order tensor minus
00310   // \return \f$ -A \f$
00311   //
00312   template<typename T>
00313   Tensor3<T>
00314   operator-(Tensor3<T> const & A)
00315   {
00316     Index const
00317     N = A.get_dimension();
00318 
00319     Tensor3<T>
00320     S(N);
00321 
00322     for (Index i = 0; i < N; ++i) {
00323       for (Index j = 0; j < N; ++j) {
00324         for (Index k = 0; k < N; ++k) {
00325           S(i,j,k) = - A(i,j,k);
00326         }
00327       }
00328     }
00329 
00330     return S;
00331   }
00332 
00333   //
00334   // 3rd-order tensor equality
00335   // Tested by components
00336   //
00337   template<typename T>
00338   inline bool
00339   operator==(Tensor3<T> const & A, Tensor3<T> const & B)
00340   {
00341     Index const
00342     N = A.get_dimension();
00343 
00344     assert(B.get_dimension() == N);
00345 
00346     for (Index i = 0; i < N; ++i) {
00347       for (Index j = 0; j < N; ++j) {
00348         for (Index k = 0; k < N; ++k) {
00349           if (A(i,j,k) != B(i,j,k)) {
00350             return false;
00351           }
00352         }
00353       }
00354     }
00355 
00356     return true;
00357   }
00358 
00359   //
00360   // 3rd-order tensor inequality
00361   // Tested by components
00362   //
00363   template<typename T>
00364   inline bool
00365   operator!=(Tensor3<T> const & A, Tensor3<T> const & B)
00366   {
00367     return !(A==B);
00368   }
00369 
00370   //
00371   // Scalar 3rd-order tensor product
00372   // \param s scalar
00373   // \param A 3rd-order tensor
00374   // \return \f$ s A \f$
00375   //
00376   template<typename T, typename S>
00377   Tensor3<T>
00378   operator*(S const & s, Tensor3<T> const & A)
00379   {
00380     Index const
00381     N = A.get_dimension();
00382 
00383     Tensor3<T>
00384     B(N);
00385 
00386     for (Index i = 0; i < N; ++i) {
00387       for (Index j = 0; j < N; ++j) {
00388         for (Index k = 0; k < N; ++k) {
00389           B(i,j,k) = s * A(i,j,k);
00390         }
00391       }
00392     }
00393 
00394     return B;
00395   }
00396 
00397   //
00398   // 3rd-order tensor scalar product
00399   // \param A 3rd-order tensor
00400   // \param s scalar
00401   // \return \f$ s A \f$
00402   //
00403   template<typename T, typename S>
00404   Tensor3<T>
00405   operator*(Tensor3<T> const & A, S const & s)
00406   {
00407     return s * A;
00408   }
00409 
00410   //
00411   // 3rd-order tensor scalar division
00412   // \param A 3rd-order tensor
00413   // \param s scalar
00414   // \return \f$ s A \f$
00415   //
00416   template<typename T, typename S>
00417   Tensor3<T>
00418   operator/(Tensor3<T> const & A, S const & s)
00419   {
00420     Index const
00421     N = A.get_dimension();
00422 
00423     Tensor3<T>
00424     B(N);
00425 
00426     for (Index i = 0; i < N; ++i) {
00427       for (Index j = 0; j < N; ++j) {
00428         for (Index k = 0; k < N; ++k) {
00429           B(i,j,k) = A(i,j,k) / s;
00430         }
00431       }
00432     }
00433 
00434     return B;
00435   }
00436 
00437   //
00438   // 3rd-order tensor vector product
00439   // \param A 3rd-order tensor
00440   // \param u vector
00441   // \return \f$ A u \f$
00442   //
00443   template<typename T>
00444   Tensor<T>
00445   dot(Tensor3<T> const & A, Vector<T> const & u)
00446   {
00447     Index const
00448     N = A.get_dimension();
00449 
00450     assert(u.get_dimension() == N);
00451 
00452     Tensor<T>
00453     B(N);
00454 
00455     for (Index j = 0; j < N; ++j) {
00456       for (Index k = 0; k < N; ++k) {
00457         T s = 0.0;
00458         for (Index i = 0; i < N; ++i) {
00459           s += A(i,j,k) * u(i);
00460         }
00461         B(j,k) = s;
00462       }
00463     }
00464 
00465     return B;
00466   }
00467 
00468   //
00469   // vector 3rd-order tensor product
00470   // \param A 3rd-order tensor
00471   // \param u vector
00472   // \return \f$ u A \f$
00473   //
00474   template<typename T>
00475   Tensor<T>
00476   dot(Vector<T> const & u, Tensor3<T> const & A)
00477   {
00478     Index const
00479     N = A.get_dimension();
00480 
00481     assert(u.get_dimension() == N);
00482 
00483     Tensor<T>
00484     B(N);
00485 
00486     for (Index i = 0; i < N; ++i) {
00487       for (Index j = 0; j < N; ++j) {
00488         T s = 0.0;
00489         for (Index k = 0; k < N; ++k) {
00490           s += A(i,j,k) * u(k);
00491         }
00492         B(i,j) = s;
00493       }
00494     }
00495 
00496     return B;
00497   }
00498 
00499 
00500   //
00501   // 3rd-order tensor vector product2 (contract 2nd index)
00502   // \param A 3rd-order tensor
00503   // \param u vector
00504   // \return \f$ A u \f$
00505   //
00506   template<typename T>
00507   Tensor<T>
00508   dot2(Tensor3<T> const & A, Vector<T> const & u)
00509   {
00510     Index const
00511     N = A.get_dimension();
00512 
00513     assert(u.get_dimension() == N);
00514 
00515     Tensor<T>
00516     B(N);
00517 
00518     for (Index i = 0; i < N; ++i) {
00519       for (Index k = 0; k < N; ++k) {
00520         T s = 0.0;
00521         for (Index j = 0; j < N; ++j) {
00522           s += A(i,j,k) * u(j);
00523         }
00524         B(i,k) = s;
00525       }
00526     }
00527 
00528     return B;
00529   }
00530 
00531   //
00532   // vector 3rd-order tensor product2 (contract 2nd index)
00533   // \param A 3rd-order tensor
00534   // \param u vector
00535   // \return \f$ u A \f$
00536   //
00537   template<typename T>
00538   Tensor<T>
00539   dot2(Vector<T> const & u, Tensor3<T> const & A)
00540   {
00541     return dot2(A, u);
00542   }
00543 
00547   template<typename T>
00548   Tensor3<T>
00549   dot(Tensor3<T> const & A, Tensor<T> const & B)
00550   {
00551     Index const
00552     N = A.get_dimension();
00553 
00554     assert(B.get_dimension() == N);
00555 
00556     Tensor3<T>
00557     C(N);
00558 
00559     for (Index i = 0; i < N; ++i) {
00560       for (Index k = 0; k < N; ++k) {
00561         for (Index j = 0; j < N; ++j) {
00562           T s = 0.0;
00563           for (Index p = 0; p < N; ++p) {
00564             s += A(i,j,p) * B(p,k);
00565           }
00566           C(i,j,k) = s;
00567         }
00568       }
00569     }
00570 
00571     return C;
00572   }
00573 
00577   template<typename T>
00578   Tensor3<T>
00579   dot(Tensor<T> const & A, Tensor3<T> const & B);
00580 
00584   template<typename T>
00585   Tensor3<T>
00586   dot2(Tensor3<T> const & A, Tensor<T> const & B);
00587 
00591   template<typename T>
00592   Tensor3<T>
00593   dot2(Tensor<T> const & A, Tensor3<T> const & B);
00594 
00595   //
00596   // 3rd-order tensor input
00597   // \param A 3rd-order tensor
00598   // \param is input stream
00599   // \return is input stream
00600   //
00601   template<typename T>
00602   std::istream &
00603   operator>>(std::istream & is, Tensor3<T> & A)
00604   {
00605     Index const
00606     N = A.get_dimension();
00607 
00608     for (Index i = 0; i < N; ++i) {
00609       for (Index j = 0; j < N; ++j) {
00610         for (Index k = 0; k < N; ++k) {
00611           is >> A(i,j,k);
00612         }
00613       }
00614     }
00615 
00616     return is;
00617   }
00618 
00619   //
00620   // 3rd-order tensor output
00621   // \param A 3rd-order tensor
00622   // \param os output stream
00623   // \return os output stream
00624   //
00625   template<typename T>
00626   std::ostream &
00627   operator<<(std::ostream & os, Tensor3<T> const & A)
00628   {
00629     Index const
00630     N = A.get_dimension();
00631 
00632     if (N == 0) {
00633       return os;
00634     }
00635 
00636     for (Index i = 0; i < N; ++i) {
00637 
00638       for (Index j = 0; j < N; ++j) {
00639 
00640         os << std::scientific << A(i,j,0);
00641 
00642         for (Index k = 1; k < N; ++k) {
00643           os << std::scientific << "," << A(i,j,k);
00644         }
00645 
00646         os << std::endl;
00647 
00648       }
00649 
00650       os << std::endl;
00651       os << std::endl;
00652 
00653     }
00654 
00655     return os;
00656   }
00657 
00658 } // namespace Intrepid
00659 
00660 #endif // Intrepid_MiniTensor_Tensor3_t_h