BidomainTissue.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 "BidomainTissue.hpp"
00030 
00031 #include "DistributedVector.hpp"
00032 #include "AxisymmetricConductivityTensors.hpp"
00033 #include "OrthotropicConductivityTensors.hpp"
00034 #include "ChastePoint.hpp"
00035 #include "AbstractChasteRegion.hpp"
00036 
00037 template <unsigned SPACE_DIM>
00038 BidomainTissue<SPACE_DIM>::BidomainTissue(
00039             AbstractCardiacCellFactory<SPACE_DIM>* pCellFactory,
00040             bool exchangeHalos)
00041     : AbstractCardiacTissue<SPACE_DIM>(pCellFactory, exchangeHalos)
00042 {
00043     CreateExtracellularConductivityTensors();
00044 }
00045 
00046 template <unsigned SPACE_DIM>
00047 BidomainTissue<SPACE_DIM>::BidomainTissue(AbstractTetrahedralMesh<SPACE_DIM,SPACE_DIM>* pMesh)
00048         :  AbstractCardiacTissue<SPACE_DIM>(pMesh)
00049 {
00050     CreateExtracellularConductivityTensors();
00051 }
00052 
00053 template <unsigned SPACE_DIM>
00054 void BidomainTissue<SPACE_DIM>::CreateExtracellularConductivityTensors()
00055 {
00056     if (this->mpConfig->IsMeshProvided() && this->mpConfig->GetLoadMesh())
00057     {
00058         assert(this->mFibreFilePathNoExtension != "");        
00059         
00060         switch (this->mpConfig->GetConductivityMedia())
00061         {
00062             case cp::media_type::Orthotropic:
00063             {
00064                 mpExtracellularConductivityTensors =  new OrthotropicConductivityTensors<SPACE_DIM,SPACE_DIM>;
00065                 FileFinder ortho_file(this->mFibreFilePathNoExtension + ".ortho", RelativeTo::AbsoluteOrCwd);
00066                 assert(ortho_file.Exists());
00067                 mpExtracellularConductivityTensors->SetFibreOrientationFile(ortho_file);
00068                 break;
00069             }
00070 
00071             case cp::media_type::Axisymmetric:
00072             {
00073                 mpExtracellularConductivityTensors =  new AxisymmetricConductivityTensors<SPACE_DIM,SPACE_DIM>;
00074                 FileFinder axi_file(this->mFibreFilePathNoExtension + ".axi", RelativeTo::AbsoluteOrCwd);
00075                 assert(axi_file.Exists());
00076                 mpExtracellularConductivityTensors->SetFibreOrientationFile(axi_file);
00077                 break;
00078             }
00079 
00080             case cp::media_type::NoFibreOrientation:
00081                 mpExtracellularConductivityTensors =  new OrthotropicConductivityTensors<SPACE_DIM,SPACE_DIM>;
00082                 break;
00083 
00084             default :
00085                 NEVER_REACHED;
00086         }
00087     }
00088     else // Slab defined in config file or SetMesh() called; no fibre orientation assumed
00089     {
00090         mpExtracellularConductivityTensors =  new OrthotropicConductivityTensors<SPACE_DIM,SPACE_DIM>;
00091     }
00092 
00093     c_vector<double, SPACE_DIM> extra_conductivities;
00094     this->mpConfig->GetExtracellularConductivities(extra_conductivities);
00095 
00096     // this definition must be here (and not inside the if statement) because SetNonConstantConductivities() will keep
00097     // a pointer to it and we don't want it to go out of scope before Init() is called
00098     unsigned num_local_elements = this->mpMesh->GetNumLocalElements();
00099     std::vector<c_vector<double, SPACE_DIM> > hetero_extra_conductivities;
00100 
00101     if (this->mpConfig->GetConductivityHeterogeneitiesProvided())
00102     {
00103         try
00104         {
00105             assert(hetero_extra_conductivities.size()==0);
00106             //initialise with the values of the default conductivity tensor 
00107             hetero_extra_conductivities.resize(num_local_elements, extra_conductivities);
00108         }
00109         catch(std::bad_alloc &r_bad_alloc)
00110         {
00111 #define COVERAGE_IGNORE
00112             std::cout << "Failed to allocate std::vector of size " << num_local_elements << std::endl;
00113             PetscTools::ReplicateException(true);
00114             throw r_bad_alloc;
00115 #undef COVERAGE_IGNORE
00116         }
00117         PetscTools::ReplicateException(false);
00118 
00119         std::vector<boost::shared_ptr<AbstractChasteRegion<SPACE_DIM> > > conductivities_heterogeneity_areas;
00120         std::vector< c_vector<double,3> > intra_h_conductivities;
00121         std::vector< c_vector<double,3> > extra_h_conductivities;
00122         HeartConfig::Instance()->GetConductivityHeterogeneities(conductivities_heterogeneity_areas,
00123                                                                 intra_h_conductivities,
00124                                                                 extra_h_conductivities);
00125 
00126         unsigned local_element_index = 0;
00127         
00128         for (typename AbstractTetrahedralMesh<SPACE_DIM,SPACE_DIM>::ElementIterator iter = (this->mpMesh)->GetElementIteratorBegin();
00129              iter != (this->mpMesh)->GetElementIteratorEnd();
00130              ++iter)
00131         {
00132             // if element centroid is contained in the region
00133             ChastePoint<SPACE_DIM> element_centroid(iter->CalculateCentroid());
00134             for (unsigned region_index=0; region_index< conductivities_heterogeneity_areas.size(); region_index++)
00135             {
00136                 // if element centroid is contained in the region
00137               if ( conductivities_heterogeneity_areas[region_index]->DoesContain( element_centroid ) )
00138                 {
00139                     //We don't use ublas vector assignment here, because we might be getting a subvector of a 3-vector
00140                     for (unsigned i=0; i<SPACE_DIM; i++)
00141                     {
00142                         hetero_extra_conductivities[local_element_index][i] = extra_h_conductivities[region_index][i];
00143                     }
00144                 }
00145             }
00146             local_element_index++;
00147         }
00148         mpExtracellularConductivityTensors->SetNonConstantConductivities(&hetero_extra_conductivities);
00149     }
00150     else
00151     {
00152         mpExtracellularConductivityTensors->SetConstantConductivities(extra_conductivities);
00153     }
00154 
00155     mpExtracellularConductivityTensors->Init(this->mpMesh);
00156 }
00157 
00158 template <unsigned SPACE_DIM>
00159 BidomainTissue<SPACE_DIM>::~BidomainTissue()
00160 {
00161     if (mpExtracellularConductivityTensors)
00162     {
00163         delete mpExtracellularConductivityTensors;
00164     }
00165 }
00166 
00167 
00168 template <unsigned SPACE_DIM>
00169 const c_matrix<double, SPACE_DIM, SPACE_DIM>& BidomainTissue<SPACE_DIM>::rGetExtracellularConductivityTensor(unsigned elementIndex)
00170 {
00171     assert(mpExtracellularConductivityTensors);
00172     if(this->mpConductivityModifier==NULL)
00173     {
00174         return (*mpExtracellularConductivityTensors)[elementIndex];
00175     }
00176     else
00177     {
00178         return this->mpConductivityModifier->rGetModifiedConductivityTensor(elementIndex, (*mpExtracellularConductivityTensors)[elementIndex]);
00179     }
00180 }
00181 
00182 
00184 // Explicit instantiation
00186 
00187 template class BidomainTissue<1>;
00188 template class BidomainTissue<2>;
00189 template class BidomainTissue<3>;
00190 
00191 
00192 // Serialization for Boost >= 1.36
00193 #include "SerializationExportWrapperForCpp.hpp"
00194 EXPORT_TEMPLATE_CLASS_SAME_DIMS(BidomainTissue)

Generated on Mon Apr 18 11:35:28 2011 for Chaste by  doxygen 1.5.5