OrthotropicConductivityTensors.cpp

00001 /*
00002 
00003 Copyright (C) University of Oxford, 2005-2009
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 SPACE_DIM>
00033 void OrthotropicConductivityTensors<SPACE_DIM>::ReadOrientationMatrixFromFile (c_matrix<double,SPACE_DIM,SPACE_DIM>& orientMatrix)
00034 {
00035     std::vector<double> items;
00036     unsigned num_elems = this->GetTokensAtNextLine(items);
00037 
00038     if (num_elems != SPACE_DIM*SPACE_DIM)
00039     {
00040         this->CloseFibreOrientationFile();
00041         EXCEPTION("Orthotropic media defined. Number of vectors in fibre orientation file and size of them should match SPACE_DIM");
00042     }
00043 
00044     for (unsigned vector_index=0; vector_index<SPACE_DIM; vector_index++)
00045     {
00046         for (unsigned component_index=0; component_index<SPACE_DIM; component_index++)
00047         {
00048             orientMatrix(component_index, vector_index) = items[vector_index*SPACE_DIM + component_index];
00049         }
00050     }
00051 }
00052 
00053 template<unsigned SPACE_DIM>
00054 void OrthotropicConductivityTensors<SPACE_DIM>::Init() throw (Exception)
00055 {
00056     if (!this->mUseNonConstantConductivities && !this->mUseFibreOrientation)
00057     {
00058         // Constant tensor for every element
00059         c_matrix<double, SPACE_DIM, SPACE_DIM> conductivity_matrix(zero_matrix<double>(SPACE_DIM,SPACE_DIM));
00060         
00061         for (unsigned dim=0; dim<SPACE_DIM; dim++)
00062         {
00063             assert(this->mConstantConductivities(dim) != DBL_MAX);
00064             conductivity_matrix(dim,dim) = this->mConstantConductivities(dim);
00065         }
00066         
00067         this->mTensors.push_back(conductivity_matrix);
00068     }
00069     else
00070     {
00071         c_matrix<double,SPACE_DIM,SPACE_DIM> orientation_matrix((identity_matrix<double>(SPACE_DIM)));
00072 
00073         if (this->mUseFibreOrientation)
00074         {
00075             this->OpenFibreOrientationFile();
00076             this->mNumElements = this->GetNumElementsFromFile();
00077         }
00078         else
00079         {
00080             this->mNumElements = this->mpNonConstantConductivities->size();
00081         }
00082 
00083         // reserve() allocates all the memory at once, more efficient than relying
00084         // on the automatic reallocation scheme.
00085         this->mTensors.reserve(this->mNumElements);
00086         
00087         c_matrix<double, SPACE_DIM, SPACE_DIM> conductivity_matrix(zero_matrix<double>(SPACE_DIM,SPACE_DIM));            
00088 
00089         for (unsigned element_index=0; element_index<this->mNumElements; element_index++)
00090         {
00091             /*
00092              *  For every element of the mesh we compute its tensor like (from
00093              * "Laminar Arrangement of VentricularMyocites Influences Electrical
00094              * Behavior of the Heart", Darren et al. 2007):
00095              *
00096              *                         [g_f  0   0 ] [a_f']
00097              *  tensor = [a_f a_l a_n] [ 0  g_l  0 ] [a_l']
00098              *                         [ 0   0  g_n] [a_n']
00099              *
00100              *              [x_i]
00101              *  where a_i = [y_i], i={f,l,n} are read from the fibre orientation file and
00102              *              [z_i]
00103              *
00104              *  g_f = fibre/longitudinal conductivity (constant or element specific)
00105              *  g_l = laminar/transverse conductivity (constant or element specific)
00106              *  g_n = normal conductivity (constant or element specific)
00107              *
00108              */
00109             if (this->mUseNonConstantConductivities)
00110             {
00111                 for (unsigned dim=0; dim<SPACE_DIM; dim++)
00112                 {
00113                     conductivity_matrix(dim,dim) = (*this->mpNonConstantConductivities)[element_index][dim];
00114                 }
00115             }
00116             else
00117             {
00118                 for (unsigned dim=0; dim<SPACE_DIM; dim++)
00119                 {
00120                     assert(this->mConstantConductivities(dim) != DBL_MAX);
00121                     conductivity_matrix(dim,dim) = this->mConstantConductivities(dim);
00122                 }                                       
00123             }
00124 
00125             if (this->mUseFibreOrientation)
00126             {
00127                 ReadOrientationMatrixFromFile(orientation_matrix);
00128             }
00129 
00130             c_matrix<double,SPACE_DIM,SPACE_DIM> temp;
00131             noalias(temp) = prod(orientation_matrix, conductivity_matrix);
00132             this->mTensors.push_back( prod(temp, trans(orientation_matrix) ) );
00133         }
00134 
00135         if (this->mUseFibreOrientation)
00136         {
00137             this->CloseFibreOrientationFile();
00138         }
00139     }
00140 
00141     this->mInitialised = true;
00142 }
00143 
00144 
00146 // Explicit instantiation
00148 
00149 template class OrthotropicConductivityTensors<1>;
00150 template class OrthotropicConductivityTensors<2>;
00151 template class OrthotropicConductivityTensors<3>;

Generated on Wed Mar 18 12:51:52 2009 for Chaste by  doxygen 1.5.5