FourthOrderTensor.hpp

00001 /*
00002 
00003 Copyright (C) University of Oxford, 2005-2010
00004 
00005 University of Oxford means the Chancellor, Masters and Scholars of the
00006 University of Oxford, having an administrative office at Wellington
00007 Square, Oxford OX1 2JD, UK.
00008 
00009 This file is part of Chaste.
00010 
00011 Chaste is free software: you can redistribute it and/or modify it
00012 under the terms of the GNU Lesser General Public License as published
00013 by the Free Software Foundation, either version 2.1 of the License, or
00014 (at your option) any later version.
00015 
00016 Chaste is distributed in the hope that it will be useful, but WITHOUT
00017 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
00018 FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
00019 License for more details. The offer of Chaste under the terms of the
00020 License is subject to the License being interpreted in accordance with
00021 English Law and subject to any action against the University of Oxford
00022 being under the jurisdiction of the English Courts.
00023 
00024 You should have received a copy of the GNU Lesser General Public License
00025 along with Chaste. If not, see <http://www.gnu.org/licenses/>.
00026 
00027 */
00028 
00029 
00030 #ifndef _FOURTHORDERTENSOR_HPP_
00031 #define _FOURTHORDERTENSOR_HPP_
00032 
00033 #include <cassert>
00034 #include <vector>
00035 
00036 #include <boost/numeric/ublas/vector.hpp>
00037 #include <boost/numeric/ublas/matrix.hpp>
00038 using namespace boost::numeric::ublas;
00039 #include "Exception.hpp"
00040 
00041 
00048 template<unsigned DIM1, unsigned DIM2, unsigned DIM3, unsigned DIM4>
00049 class FourthOrderTensor
00050 {
00051 private:
00052     std::vector<double> mData;  
00060     unsigned GetVectorIndex(unsigned M, unsigned N, unsigned P, unsigned Q)
00061     {
00062         assert(M<DIM1);
00063         assert(N<DIM2);
00064         assert(P<DIM3);
00065         assert(Q<DIM4);
00066         return M + DIM1*N + DIM1*DIM2*P + DIM1*DIM2*DIM3*Q;
00067     }
00068 
00069 public:
00070 
00074     FourthOrderTensor();
00075 
00084     template<unsigned CONTRACTED_DIM>
00085     void SetAsContractionOnFirstDimension(const c_matrix<double,DIM1,CONTRACTED_DIM>& rMatrix, FourthOrderTensor<CONTRACTED_DIM,DIM2,DIM3,DIM4>& rTensor);
00086 
00087 
00096     template<unsigned CONTRACTED_DIM>
00097     void SetAsContractionOnSecondDimension(const c_matrix<double,DIM2,CONTRACTED_DIM>& rMatrix, FourthOrderTensor<DIM1,CONTRACTED_DIM,DIM3,DIM4>& rTensor);
00098 
00107     template<unsigned CONTRACTED_DIM>
00108     void SetAsContractionOnThirdDimension(const c_matrix<double,DIM3,CONTRACTED_DIM>& rMatrix, FourthOrderTensor<DIM1,DIM2,CONTRACTED_DIM,DIM4>& rTensor);
00109 
00119     template<unsigned CONTRACTED_DIM>
00120     void SetAsContractionOnFourthDimension(const c_matrix<double,DIM4,CONTRACTED_DIM>& rMatrix, FourthOrderTensor<DIM1,DIM2,DIM3,CONTRACTED_DIM>& rTensor);
00121 
00122 
00131     double& operator()(unsigned M, unsigned N, unsigned P, unsigned Q);
00132 
00136     void Zero();
00137     
00141     std::vector<double>& rGetData()
00142     {
00143         return mData;
00144     }
00145 };
00146 
00147 
00148 
00149 
00151 // Implementation (lots of possibilities for the dimensions so no point with expl instantiation)
00153 
00154 
00155 template<unsigned DIM1, unsigned DIM2, unsigned DIM3, unsigned DIM4>
00156 FourthOrderTensor<DIM1,DIM2,DIM3,DIM4>::FourthOrderTensor()
00157 {
00158     unsigned size = DIM1*DIM2*DIM3*DIM4;
00159     mData.resize(size, 0.0);
00160 }
00161 
00162 template<unsigned DIM1, unsigned DIM2, unsigned DIM3, unsigned DIM4>
00163 template<unsigned CONTRACTED_DIM>
00164 void FourthOrderTensor<DIM1,DIM2,DIM3,DIM4>::SetAsContractionOnFirstDimension(const c_matrix<double,DIM1,CONTRACTED_DIM>& rMatrix, FourthOrderTensor<CONTRACTED_DIM,DIM2,DIM3,DIM4>& rTensor)
00165 {
00166     Zero();
00167 
00168     std::vector<double>::iterator iter = mData.begin();
00169     std::vector<double>::iterator other_tensor_iter = rTensor.rGetData().begin();
00170 
00171     for(unsigned d=0; d<DIM4; d++)
00172     {
00173         for(unsigned c=0; c<DIM3; c++)
00174         {
00175             for(unsigned b=0; b<DIM2; b++)
00176             {
00177                 for(unsigned a=0; a<DIM1; a++)
00178                 {
00179                     for(unsigned N=0; N<CONTRACTED_DIM; N++)
00180                     {
00181                         // The following just does
00182                         //
00183                         // mData[GetVectorIndex(a,b,c,d)] += rMatrix(a,N) * rTensor(N,b,c,d);
00184                         //
00185                         // but more efficiently using iterators into the data vector, not
00186                         // using random access
00187                         *iter += rMatrix(a,N) * *other_tensor_iter;
00188                         other_tensor_iter++;
00189                     }
00190                     
00191                     iter++;
00192                     
00193                     if(a != DIM1-1)
00194                     {
00195                         other_tensor_iter -= CONTRACTED_DIM;
00196                     }
00197                 }
00198             }
00199         }
00200     }
00201 }
00202  
00203 
00204 template<unsigned DIM1, unsigned DIM2, unsigned DIM3, unsigned DIM4>
00205 template<unsigned CONTRACTED_DIM>
00206 void FourthOrderTensor<DIM1,DIM2,DIM3,DIM4>::SetAsContractionOnSecondDimension(const c_matrix<double,DIM2,CONTRACTED_DIM>& rMatrix, FourthOrderTensor<DIM1,CONTRACTED_DIM,DIM3,DIM4>& rTensor)
00207 {
00208     Zero();
00209 
00210     std::vector<double>::iterator iter = mData.begin();
00211     std::vector<double>::iterator other_tensor_iter = rTensor.rGetData().begin();
00212 
00213     for(unsigned d=0; d<DIM4; d++)
00214     {
00215         for(unsigned c=0; c<DIM3; c++)
00216         {
00217             for(unsigned b=0; b<DIM2; b++)
00218             {
00219                 for(unsigned N=0; N<CONTRACTED_DIM; N++)
00220                 {
00221                     for(unsigned a=0; a<DIM1; a++)
00222                     {
00223                         // The following just does
00224                         //
00225                         // mData[GetVectorIndex(a,b,c,d)] += rMatrix(b,N) * rTensor(a,N,c,d);
00226                         //
00227                         // but more efficiently using iterators into the data vector, not
00228                         // using random access
00229                         *iter += rMatrix(b,N) * *other_tensor_iter;
00230                         iter++;
00231                         other_tensor_iter++;
00232                     }
00233                     
00234                     if(N != CONTRACTED_DIM-1)
00235                     {
00236                         iter -= DIM1;
00237                     }
00238                 }
00239                 if(b != DIM2-1)
00240                 {
00241                     other_tensor_iter -= CONTRACTED_DIM*DIM1;
00242                 }
00243             }
00244         }
00245     }
00246 }
00247 
00248 
00249 template<unsigned DIM1, unsigned DIM2, unsigned DIM3, unsigned DIM4>
00250 template<unsigned CONTRACTED_DIM>
00251 void FourthOrderTensor<DIM1,DIM2,DIM3,DIM4>::SetAsContractionOnThirdDimension(const c_matrix<double,DIM3,CONTRACTED_DIM>& rMatrix, FourthOrderTensor<DIM1,DIM2,CONTRACTED_DIM,DIM4>& rTensor)
00252 {
00253     Zero();
00254     
00255     std::vector<double>::iterator iter = mData.begin();
00256     std::vector<double>::iterator other_tensor_iter = rTensor.rGetData().begin();
00257     
00258     for(unsigned d=0; d<DIM4; d++)
00259     {
00260         for(unsigned c=0; c<DIM3; c++)
00261         {
00262             for(unsigned N=0; N<CONTRACTED_DIM; N++)
00263             {
00264                 for(unsigned b=0; b<DIM2; b++)
00265                 {
00266                     for(unsigned a=0; a<DIM1; a++)
00267                     {
00268                         // The following just does
00269                         //
00270                         // mData[GetVectorIndex(a,b,c,d)] += rMatrix(c,N) * rTensor(a,b,N,d);
00271                         //
00272                         // but more efficiently using iterators into the data vector, not
00273                         // using random access
00274                         *iter += rMatrix(c,N) * *other_tensor_iter;
00275                         iter++;
00276                         other_tensor_iter++;
00277                     }
00278                 }
00279 
00280                 if(N != CONTRACTED_DIM-1)
00281                 {
00282                     iter -= DIM1*DIM2;
00283                 }             
00284             }
00285 
00286             if(c != DIM3-1)
00287             {
00288                 other_tensor_iter -= CONTRACTED_DIM*DIM1*DIM2;
00289             }
00290         }
00291     }
00292 }
00293 
00294 
00295 template<unsigned DIM1, unsigned DIM2, unsigned DIM3, unsigned DIM4>
00296 template<unsigned CONTRACTED_DIM>
00297 void FourthOrderTensor<DIM1,DIM2,DIM3,DIM4>::SetAsContractionOnFourthDimension(const c_matrix<double,DIM4,CONTRACTED_DIM>& rMatrix, FourthOrderTensor<DIM1,DIM2,DIM3,CONTRACTED_DIM>& rTensor)
00298 {
00299     Zero();
00300    
00301     std::vector<double>::iterator iter = mData.begin();
00302     std::vector<double>::iterator other_tensor_iter = rTensor.rGetData().begin();
00303 
00304     for(unsigned d=0; d<DIM4; d++)
00305     {
00306         for(unsigned N=0; N<CONTRACTED_DIM; N++)
00307         {
00308             for(unsigned c=0; c<DIM3; c++)
00309             {
00310                 for(unsigned b=0; b<DIM2; b++)
00311                 {
00312                     for(unsigned a=0; a<DIM1; a++)
00313                     {
00314                         // The following just does
00315                         //
00316                         // mData[GetVectorIndex(a,b,c,d)] += rMatrix(d,N) * rTensor(a,b,c,N);
00317                         //
00318                         // but more efficiently using iterators into the data vector, not
00319                         // using random access
00320                         *iter += rMatrix(d,N) * *other_tensor_iter;
00321                         
00322                         iter++;
00323                         other_tensor_iter++;
00324                     }
00325                 }
00326             }
00327     
00328             if(N != CONTRACTED_DIM-1)
00329             {
00330                 iter-= DIM1*DIM2*DIM3;
00331             }  
00332         }
00333 
00334         other_tensor_iter -= CONTRACTED_DIM*DIM1*DIM2*DIM3;
00335     }
00336 }
00337 
00338 template<unsigned DIM1, unsigned DIM2, unsigned DIM3, unsigned DIM4>
00339 double& FourthOrderTensor<DIM1,DIM2,DIM3,DIM4>::operator()(unsigned M, unsigned N, unsigned P, unsigned Q)
00340 {
00341     assert(M<DIM1);
00342     assert(N<DIM2);
00343     assert(P<DIM3);
00344     assert(Q<DIM4);
00345 
00346     return mData[GetVectorIndex(M,N,P,Q)];
00347 }
00348 
00349 template<unsigned DIM1, unsigned DIM2, unsigned DIM3, unsigned DIM4>
00350 void FourthOrderTensor<DIM1,DIM2,DIM3,DIM4>::Zero()
00351 {
00352     for (unsigned i=0; i<mData.size(); i++)
00353     {
00354         mData[i] = 0.0;
00355     }
00356 }
00357 
00358 #endif //_FOURTHORDERTENSOR_HPP_

Generated on Mon Nov 1 12:35:17 2010 for Chaste by  doxygen 1.5.5