OrthotropicConductivityTensors.cpp

00001 /*
00002 
00003 Copyright (C) University of Oxford, 2005-2011
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 #include "OrthotropicConductivityTensors.hpp"
00030 #include "Exception.hpp"
00031 
00032 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00033 void OrthotropicConductivityTensors<ELEMENT_DIM, SPACE_DIM>::Init(AbstractTetrahedralMesh<ELEMENT_DIM, SPACE_DIM> *pMesh) throw (Exception)
00034 {
00035     this->mpMesh = pMesh;
00036 
00037     if (!this->mUseNonConstantConductivities && !this->mUseFibreOrientation)
00038     {
00039         // Constant tensor for every element
00040         c_matrix<double, SPACE_DIM, SPACE_DIM> conductivity_matrix(zero_matrix<double>(SPACE_DIM,SPACE_DIM));
00041 
00042         for (unsigned dim=0; dim<SPACE_DIM; dim++)
00043         {
00044             assert(this->mConstantConductivities(dim) != DBL_MAX);
00045             conductivity_matrix(dim,dim) = this->mConstantConductivities(dim);
00046         }
00047 
00048         this->mTensors.push_back(conductivity_matrix);
00049     }
00050     else
00051     {
00052         c_matrix<double,SPACE_DIM,SPACE_DIM> orientation_matrix((identity_matrix<double>(SPACE_DIM)));
00053 
00054         if (this->mUseFibreOrientation)
00055         {
00056             // open file
00057             this->mFileReader.reset(new FibreReader<SPACE_DIM>(this->mFibreOrientationFile, ORTHO));
00058             if(this->mFileReader->GetNumLinesOfData() != this->mpMesh->GetNumElements())
00059             {
00060                 EXCEPTION("The size of the fibre file does not match the number of elements in the mesh");
00061             }
00062         }
00063  
00064         if (this->mUseNonConstantConductivities)
00065         {
00066             if(this->mpNonConstantConductivities->size() != this->mpMesh->GetNumLocalElements())
00067             {
00068                 EXCEPTION("The size of the conductivities vector does not match the number of elements in the mesh");
00069             }
00070         }
00071 
00072         // reserve() allocates all the memory at once, more efficient than relying
00073         // on the automatic reallocation scheme.
00074         this->mTensors.reserve(this->mpMesh->GetNumLocalElements());
00075 
00076         c_matrix<double, SPACE_DIM, SPACE_DIM> conductivity_matrix(zero_matrix<double>(SPACE_DIM,SPACE_DIM));
00077 
00078         unsigned local_element_index = 0;
00079 
00080         int previous_global_index = -1;
00081 
00082         for (typename AbstractTetrahedralMesh<ELEMENT_DIM,SPACE_DIM>::ElementIterator it = this->mpMesh->GetElementIteratorBegin();
00083              it != this->mpMesh->GetElementIteratorEnd();
00084              ++it)
00085         {
00086             if (this->mUseFibreOrientation)
00087             {
00088                 int current_fibre_global_index = (int) it->GetIndex();
00089 
00090                 // Assumption about ElementIterator returning elements in ascending order is wrong
00091                 // if this fails
00092                 assert(current_fibre_global_index > previous_global_index);
00093 
00094                 
00095                 for (int fibre_index=previous_global_index; fibre_index<current_fibre_global_index-1; fibre_index++)
00096                 {
00097                     this->mFileReader->GetNextFibreSheetAndNormalMatrix(orientation_matrix);            
00098                 }
00099                 
00100                 previous_global_index = current_fibre_global_index;
00101             }
00102 
00103             /*
00104              *  For every element of the mesh we compute its tensor like (from
00105              * "Laminar Arrangement of VentricularMyocites Influences Electrical
00106              * Behavior of the Heart", Darren et al. 2007):
00107              *
00108              *                         [g_f  0   0 ] [a_f']
00109              *  tensor = [a_f a_l a_n] [ 0  g_l  0 ] [a_l']
00110              *                         [ 0   0  g_n] [a_n']
00111              *
00112              *              [x_i]
00113              *  where a_i = [y_i], i={f,l,n} are read from the fibre orientation file and
00114              *              [z_i]
00115              *
00116              *  g_f = fibre/longitudinal conductivity (constant or element specific)
00117              *  g_l = laminar/transverse conductivity (constant or element specific)
00118              *  g_n = normal conductivity (constant or element specific)
00119              *
00120              */
00121             if (this->mUseNonConstantConductivities)
00122             {
00123                 for (unsigned dim=0; dim<SPACE_DIM; dim++)
00124                 {
00125                     conductivity_matrix(dim,dim) = (*this->mpNonConstantConductivities)[local_element_index][dim];
00126                 }
00127             }
00128             else
00129             {
00130                 for (unsigned dim=0; dim<SPACE_DIM; dim++)
00131                 {
00132                     assert(this->mConstantConductivities(dim) != DBL_MAX);
00133                     conductivity_matrix(dim,dim) = this->mConstantConductivities(dim);
00134                 }
00135             }
00136 
00137             if (this->mUseFibreOrientation)
00138             {
00139                 this->mFileReader->GetNextFibreSheetAndNormalMatrix(orientation_matrix);
00140             }
00141 
00142             c_matrix<double,SPACE_DIM,SPACE_DIM> temp;
00143             noalias(temp) = prod(orientation_matrix, conductivity_matrix);
00144             this->mTensors.push_back( prod(temp, trans(orientation_matrix) ) );
00145             
00146             local_element_index++;
00147         }
00148         assert(this->mTensors.size() == this->mpMesh->GetNumLocalElements());
00149         assert(this->mTensors.size() == local_element_index);
00150 
00151         if (this->mUseFibreOrientation)
00152         {
00153             // close fibre file
00154             this->mFileReader.reset();
00155         }
00156     }
00157 
00158     this->mInitialised = true;
00159 }
00160 
00161 
00163 // Explicit instantiation
00165 
00166 template class OrthotropicConductivityTensors<1,1>;
00167 template class OrthotropicConductivityTensors<1,2>;
00168 template class OrthotropicConductivityTensors<1,3>;
00169 template class OrthotropicConductivityTensors<2,2>;
00170 template class OrthotropicConductivityTensors<2,3>;
00171 template class OrthotropicConductivityTensors<3,3>;

Generated on Tue May 31 14:31:41 2011 for Chaste by  doxygen 1.5.5