BidomainPde.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 "BidomainPde.hpp"
00030 
00031 #include "DistributedVector.hpp"
00032 #include "AxisymmetricConductivityTensors.hpp"
00033 #include "OrthotropicConductivityTensors.hpp"
00034 #include "ChastePoint.hpp"
00035 #include "ChasteCuboid.hpp"
00036 
00037 
00038 template <unsigned SPACE_DIM>
00039 BidomainPde<SPACE_DIM>::BidomainPde(
00040             AbstractCardiacCellFactory<SPACE_DIM>* pCellFactory)
00041     : AbstractCardiacPde<SPACE_DIM>(pCellFactory, 2 /*mStride*/)
00042 {
00043     mExtracellularStimulusCacheReplicated.resize( pCellFactory->GetNumberOfCells() );
00044 
00045     if (this->mpConfig->GetIsMeshProvided() && this->mpConfig->GetLoadMesh())
00046     {
00047         switch (this->mpConfig->GetConductivityMedia())
00048         {
00049             case media_type::Orthotropic:
00050                 mpExtracellularConductivityTensors =  new OrthotropicConductivityTensors<SPACE_DIM>;
00051                 mpExtracellularConductivityTensors->SetFibreOrientationFile(this->mpConfig->GetMeshName() + ".ortho");                
00052                 break;
00053                 
00054             case media_type::Axisymmetric:
00055                 mpExtracellularConductivityTensors =  new AxisymmetricConductivityTensors<SPACE_DIM>;
00056                 mpExtracellularConductivityTensors->SetFibreOrientationFile(this->mpConfig->GetMeshName() + ".axi");
00057                 break;
00058 
00059             case media_type::NoFibreOrientation:
00060                 mpExtracellularConductivityTensors =  new OrthotropicConductivityTensors<SPACE_DIM>;                
00061                 break;
00062                 
00063             default :
00064                 NEVER_REACHED;                                
00065         }            
00066     }
00067     else // Slab defined in config file or SetMesh() called; no fibre orientation assumed
00068     {
00069         mpExtracellularConductivityTensors =  new OrthotropicConductivityTensors<SPACE_DIM>;           
00070     }
00071 
00072     c_vector<double, SPACE_DIM> extra_conductivities;
00073     this->mpConfig->GetExtracellularConductivities(extra_conductivities);    
00074 
00075     // this definition must be here (and not inside the if statement) because SetNonConstantConductivities() will keep 
00076     // a pointer to it and we don't want it to go out of scope before Init() is called
00077     unsigned num_elements = pCellFactory->GetMesh()->GetNumElements();                        
00078     std::vector<c_vector<double, SPACE_DIM> > hetero_extra_conductivities(num_elements); 
00079 
00080     if (this->mpConfig->GetConductivityHeterogeneitiesProvided())
00081     {
00082         std::vector<ChasteCuboid> conductivities_heterogeneity_areas;
00083         std::vector< c_vector<double,3> > intra_h_conductivities;
00084         std::vector< c_vector<double,3> > extra_h_conductivities;
00085         HeartConfig::Instance()->GetConductivityHeterogeneities(conductivities_heterogeneity_areas,
00086                                                                 intra_h_conductivities,
00087                                                                 extra_h_conductivities);
00088                     
00089         for (unsigned element_index=0; element_index<num_elements; element_index++)
00090         {
00091             for (unsigned region_index=0; region_index< conductivities_heterogeneity_areas.size(); region_index++)
00092             {
00093                 // if element centroid is contained in the region
00094                 ChastePoint<SPACE_DIM> element_centroid(pCellFactory->GetMesh()->GetElement(element_index)->CalculateCentroid());
00095                 if ( conductivities_heterogeneity_areas[region_index].DoesContain( element_centroid ) )
00096                 {
00097                     hetero_extra_conductivities[element_index] = extra_h_conductivities[region_index];
00098                 }
00099                 else
00100                 {
00101                     hetero_extra_conductivities[element_index] = extra_conductivities;                        
00102                 }
00103             }
00104         }
00105         
00106         mpExtracellularConductivityTensors->SetNonConstantConductivities(&hetero_extra_conductivities);            
00107     }
00108     else
00109     {                       
00110         mpExtracellularConductivityTensors->SetConstantConductivities(extra_conductivities);
00111     }
00112     
00113     mpExtracellularConductivityTensors->Init();
00114 }
00115 
00116 template <unsigned SPACE_DIM>
00117 BidomainPde<SPACE_DIM>::~BidomainPde()
00118 {
00119     delete mpExtracellularConductivityTensors;
00120 }
00121 
00122 
00123 template <unsigned SPACE_DIM>
00124 const c_matrix<double, SPACE_DIM, SPACE_DIM>& BidomainPde<SPACE_DIM>::rGetExtracellularConductivityTensor(unsigned elementIndex)
00125 {
00126     assert(mpExtracellularConductivityTensors);
00127     return (*mpExtracellularConductivityTensors)[elementIndex];
00128 }
00129 
00130 
00132 // Explicit instantiation
00134 
00135 template class BidomainPde<1>;
00136 template class BidomainPde<2>;
00137 template class BidomainPde<3>;

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