AbstractCardiacTissue.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 "AbstractCardiacTissue.hpp"
00030 
00031 #include "DistributedVector.hpp"
00032 #include "AxisymmetricConductivityTensors.hpp"
00033 #include "OrthotropicConductivityTensors.hpp"
00034 #include "ChastePoint.hpp"
00035 #include "AbstractChasteRegion.hpp"
00036 #include "HeartEventHandler.hpp"
00037 #include "PetscTools.hpp"
00038 #include "PetscVecTools.hpp"
00039 
00040 template <unsigned ELEMENT_DIM,unsigned SPACE_DIM>
00041 AbstractCardiacTissue<ELEMENT_DIM,SPACE_DIM>::AbstractCardiacTissue(
00042             AbstractCardiacCellFactory<ELEMENT_DIM,SPACE_DIM>* pCellFactory,
00043             bool exchangeHalos)
00044     : mpMesh(pCellFactory->GetMesh()),
00045       mDoCacheReplication(true),
00046       mpDistributedVectorFactory(mpMesh->GetDistributedVectorFactory()),
00047       mMeshUnarchived(false),
00048       mpConductivityModifier(NULL),
00049       mExchangeHalos(exchangeHalos)
00050 {
00051     //This constructor is called from the Initialise() method of the CardiacProblem class
00052     assert(pCellFactory != NULL);
00053     assert(pCellFactory->GetMesh() != NULL);
00054 
00055     if (PetscTools::IsSequential())
00056     {
00057         //Remove the request for a halo exchange
00058         mExchangeHalos = false;
00059     }
00060 
00061     unsigned num_local_nodes = mpDistributedVectorFactory->GetLocalOwnership();
00062     unsigned ownership_range_low = mpDistributedVectorFactory->GetLow();
00063     mCellsDistributed.resize(num_local_nodes);
00064 
00065     try
00066     {
00067         for (unsigned local_index = 0; local_index < num_local_nodes; local_index++)
00068         {
00069             unsigned global_index = ownership_range_low + local_index;
00070             mCellsDistributed[local_index] = pCellFactory->CreateCardiacCellForNode(global_index);
00071             mCellsDistributed[local_index]->SetUsedInTissueSimulation();
00072         }
00073 
00074         pCellFactory->FinaliseCellCreation(&mCellsDistributed,
00075                                            mpDistributedVectorFactory->GetLow(),
00076                                            mpDistributedVectorFactory->GetHigh());
00077     }
00078     catch (const Exception& e)
00079     {
00080         // This catch statement is quite tricky to cover, but it is actually done in TestCardiacSimulation::TestMono1dSodiumBlockBySettingNamedParameter
00081 
00082         // Errors thrown creating cells will often be process-specific
00083         PetscTools::ReplicateException(true);
00084 
00085         // Delete cells
00086         // Should really do this for other processes too, but this is all we need
00087         // to get memory testing to pass, and leaking when we're about to die isn't
00088         // that bad!
00089         for (std::vector<AbstractCardiacCell*>::iterator cell_iterator = mCellsDistributed.begin();
00090              cell_iterator != mCellsDistributed.end();
00091              ++cell_iterator)
00092         {
00093             delete (*cell_iterator);
00094         }
00095 
00096         throw e;
00097     }
00098     PetscTools::ReplicateException(false);
00099 
00100     // Halo nodes (if required)
00101     SetUpHaloCells(pCellFactory);
00102 
00103     HeartEventHandler::BeginEvent(HeartEventHandler::COMMUNICATION);
00104     mIionicCacheReplicated.Resize( pCellFactory->GetNumberOfCells() );
00105     mIntracellularStimulusCacheReplicated.Resize( pCellFactory->GetNumberOfCells() );
00106     HeartEventHandler::EndEvent(HeartEventHandler::COMMUNICATION);
00107 
00108     if(HeartConfig::Instance()->IsMeshProvided() && HeartConfig::Instance()->GetLoadMesh())
00109     {
00110         mFibreFilePathNoExtension = "./" + HeartConfig::Instance()->GetMeshName();
00111     }
00112     else
00113     {
00114         // As of 10671 fibre orientation can only be defined when loading a mesh from disc.
00115         mFibreFilePathNoExtension = "";
00116     }
00117     CreateIntracellularConductivityTensor();
00118 }
00119 
00120 // Constructor used for archiving
00121 template <unsigned ELEMENT_DIM,unsigned SPACE_DIM>
00122 AbstractCardiacTissue<ELEMENT_DIM,SPACE_DIM>::AbstractCardiacTissue(AbstractTetrahedralMesh<ELEMENT_DIM,SPACE_DIM>* pMesh)
00123     : mpMesh(pMesh),
00124       mDoCacheReplication(true),
00125       mpDistributedVectorFactory(mpMesh->GetDistributedVectorFactory()),
00126       mMeshUnarchived(true),
00127       mExchangeHalos(false)
00128 {
00129     mIionicCacheReplicated.Resize(mpDistributedVectorFactory->GetProblemSize());
00130     mIntracellularStimulusCacheReplicated.Resize(mpDistributedVectorFactory->GetProblemSize());
00131 
00132     mFibreFilePathNoExtension = ArchiveLocationInfo::GetArchiveDirectory() + ArchiveLocationInfo::GetMeshFilename();
00133     CreateIntracellularConductivityTensor();
00134 }
00135 
00136 template <unsigned ELEMENT_DIM,unsigned SPACE_DIM>
00137 AbstractCardiacTissue<ELEMENT_DIM,SPACE_DIM>::~AbstractCardiacTissue()
00138 {
00139     // Delete cells
00140     for (std::vector<AbstractCardiacCell*>::iterator cell_iterator = mCellsDistributed.begin();
00141          cell_iterator != mCellsDistributed.end();
00142          ++cell_iterator)
00143     {
00144         delete (*cell_iterator);
00145     }
00146 
00147     // Delete cells for halo nodes
00148     for (std::vector<AbstractCardiacCell*>::iterator cell_iterator = mHaloCellsDistributed.begin();
00149          cell_iterator != mHaloCellsDistributed.end();
00150          ++cell_iterator)
00151     {
00152         delete (*cell_iterator);
00153     }
00154 
00155     delete mpIntracellularConductivityTensors;
00156 
00157     // If the mesh was unarchived we need to free it explicitly.
00158     if (mMeshUnarchived)
00159     {
00160         delete mpMesh;
00161     }
00162 }
00163 
00164 template <unsigned ELEMENT_DIM,unsigned SPACE_DIM>
00165 void AbstractCardiacTissue<ELEMENT_DIM,SPACE_DIM>::CreateIntracellularConductivityTensor()
00166 {
00167     HeartEventHandler::BeginEvent(HeartEventHandler::READ_MESH);
00168     mpConfig = HeartConfig::Instance();
00169 
00170     if (mpConfig->IsMeshProvided() && mpConfig->GetLoadMesh())
00171     {
00172         assert(mFibreFilePathNoExtension != "");
00173 
00174         switch (mpConfig->GetConductivityMedia())
00175         {
00176             case cp::media_type::Orthotropic:
00177             {
00178                 mpIntracellularConductivityTensors = new OrthotropicConductivityTensors<ELEMENT_DIM,SPACE_DIM>;
00179                 FileFinder ortho_file(mFibreFilePathNoExtension + ".ortho", RelativeTo::AbsoluteOrCwd);
00180                 assert(ortho_file.Exists());
00181                 mpIntracellularConductivityTensors->SetFibreOrientationFile(ortho_file);
00182                 break;
00183             }
00184 
00185             case cp::media_type::Axisymmetric:
00186             {
00187                 mpIntracellularConductivityTensors = new AxisymmetricConductivityTensors<ELEMENT_DIM,SPACE_DIM>;
00188                 FileFinder axi_file(mFibreFilePathNoExtension + ".axi", RelativeTo::AbsoluteOrCwd);
00189                 assert(axi_file.Exists());
00190                 mpIntracellularConductivityTensors->SetFibreOrientationFile(axi_file);
00191                 break;
00192             }
00193 
00194             case cp::media_type::NoFibreOrientation:
00196                 mpIntracellularConductivityTensors = new OrthotropicConductivityTensors<ELEMENT_DIM,SPACE_DIM>;
00197                 break;
00198 
00199             default :
00200                 NEVER_REACHED;
00201         }
00202     }
00203     else // Slab defined in config file or SetMesh() called; no fibre orientation assumed
00204     {
00206         mpIntracellularConductivityTensors = new OrthotropicConductivityTensors<ELEMENT_DIM,SPACE_DIM>;
00207     }
00208 
00209     c_vector<double, SPACE_DIM> intra_conductivities;
00210     mpConfig->GetIntracellularConductivities(intra_conductivities);
00211 
00212     // this definition must be here (and not inside the if statement) because SetNonConstantConductivities() will keep
00213     // a pointer to it and we don't want it to go out of scope before Init() is called
00214     unsigned num_local_elements = mpMesh->GetNumLocalElements();
00215     std::vector<c_vector<double, SPACE_DIM> > hetero_intra_conductivities;
00216 
00217     if (mpConfig->GetConductivityHeterogeneitiesProvided())
00218     {
00219         try
00220         {
00221             assert(hetero_intra_conductivities.size()==0);
00222             hetero_intra_conductivities.resize(num_local_elements, intra_conductivities);
00223         }
00224         catch(std::bad_alloc &r_bad_alloc)
00225         {
00226 #define COVERAGE_IGNORE
00227             std::cout << "Failed to allocate std::vector of size " << num_local_elements << std::endl;
00228             PetscTools::ReplicateException(true);
00229             throw r_bad_alloc;
00230 #undef COVERAGE_IGNORE
00231         }
00232         PetscTools::ReplicateException(false);
00233 
00234         std::vector<boost::shared_ptr<AbstractChasteRegion<SPACE_DIM> > > conductivities_heterogeneity_areas;
00235         std::vector< c_vector<double,3> > intra_h_conductivities;
00236         std::vector< c_vector<double,3> > extra_h_conductivities;
00237         HeartConfig::Instance()->GetConductivityHeterogeneities(conductivities_heterogeneity_areas,
00238                                                                 intra_h_conductivities,
00239                                                                 extra_h_conductivities);
00240 
00241         unsigned local_element_index = 0;
00242 
00243         for (typename AbstractTetrahedralMesh<ELEMENT_DIM,SPACE_DIM>::ElementIterator it = mpMesh->GetElementIteratorBegin();
00244              it != mpMesh->GetElementIteratorEnd();
00245              ++it)
00246         {
00247 //            unsigned element_index = it->GetIndex();
00248             // if element centroid is contained in the region
00249             ChastePoint<SPACE_DIM> element_centroid(it->CalculateCentroid());
00250             for (unsigned region_index=0; region_index< conductivities_heterogeneity_areas.size(); region_index++)
00251             {
00252                 if ( conductivities_heterogeneity_areas[region_index]->DoesContain(element_centroid) )
00253                 {
00254                     //We don't use ublas vector assignment here, because we might be getting a subvector of a 3-vector
00255                     for (unsigned i=0; i<SPACE_DIM; i++)
00256                     {
00257                         hetero_intra_conductivities[local_element_index][i] = intra_h_conductivities[region_index][i];
00258                     }
00259                 }
00260             }
00261             local_element_index++;
00262         }
00263 
00264         mpIntracellularConductivityTensors->SetNonConstantConductivities(&hetero_intra_conductivities);
00265     }
00266     else
00267     {
00268         mpIntracellularConductivityTensors->SetConstantConductivities(intra_conductivities);
00269     }
00270 
00271     mpIntracellularConductivityTensors->Init(this->mpMesh);
00272     HeartEventHandler::EndEvent(HeartEventHandler::READ_MESH);
00273 }
00274 
00275 
00276 template <unsigned ELEMENT_DIM,unsigned SPACE_DIM>
00277 void AbstractCardiacTissue<ELEMENT_DIM,SPACE_DIM>::SetCacheReplication(bool doCacheReplication)
00278 {
00279     mDoCacheReplication = doCacheReplication;
00280 }
00281 
00282 template <unsigned ELEMENT_DIM,unsigned SPACE_DIM>
00283 bool AbstractCardiacTissue<ELEMENT_DIM,SPACE_DIM>::GetDoCacheReplication()
00284 {
00285     return mDoCacheReplication;
00286 }
00287 
00288 template <unsigned ELEMENT_DIM,unsigned SPACE_DIM>
00289 const c_matrix<double, SPACE_DIM, SPACE_DIM>& AbstractCardiacTissue<ELEMENT_DIM,SPACE_DIM>::rGetIntracellularConductivityTensor(unsigned elementIndex)
00290 {
00291     assert( mpIntracellularConductivityTensors);
00292     if(mpConductivityModifier==NULL)
00293     {
00294         return (*mpIntracellularConductivityTensors)[elementIndex];
00295     }
00296     else
00297     {
00298         return mpConductivityModifier->rGetModifiedConductivityTensor(elementIndex, (*mpIntracellularConductivityTensors)[elementIndex]);
00299     }
00300 }
00301 
00302 template <unsigned ELEMENT_DIM,unsigned SPACE_DIM>
00303 const c_matrix<double, SPACE_DIM, SPACE_DIM>& AbstractCardiacTissue<ELEMENT_DIM,SPACE_DIM>::rGetExtracellularConductivityTensor(unsigned elementIndex)
00304 {
00305      EXCEPTION("Monodomain tissues do not have extracellular conductivity tensors.");
00306 }
00307 
00308 template <unsigned ELEMENT_DIM,unsigned SPACE_DIM>
00309 AbstractCardiacCell* AbstractCardiacTissue<ELEMENT_DIM,SPACE_DIM>::GetCardiacCell( unsigned globalIndex )
00310 {
00311     assert(mpDistributedVectorFactory->GetLow() <= globalIndex &&
00312            globalIndex < mpDistributedVectorFactory->GetHigh());
00313     return mCellsDistributed[globalIndex - mpDistributedVectorFactory->GetLow()];
00314 }
00315 
00316 template <unsigned ELEMENT_DIM,unsigned SPACE_DIM>
00317 AbstractCardiacCell* AbstractCardiacTissue<ELEMENT_DIM,SPACE_DIM>::GetCardiacCellOrHaloCell( unsigned globalIndex )
00318 {
00319     std::map<unsigned, unsigned>::const_iterator node_position;
00320     // First search the halo
00321     if ((node_position=mHaloGlobalToLocalIndexMap.find(globalIndex)) != mHaloGlobalToLocalIndexMap.end())
00322     {
00323         //Found a halo node
00324         return mHaloCellsDistributed[node_position->second];
00325     }
00326     // Then search the owned node
00327     if ( mpDistributedVectorFactory->IsGlobalIndexLocal(globalIndex)  )
00328     {
00329         //Found an owned node
00330         return mCellsDistributed[globalIndex - mpDistributedVectorFactory->GetLow()];
00331     }
00332     // Not here
00333     std::stringstream message;
00334     message << "Requested node/halo " << globalIndex << " does not belong to processor " << PetscTools::GetMyRank();
00335     EXCEPTION(message.str().c_str());
00336 }
00337 
00338 
00339 template <unsigned ELEMENT_DIM,unsigned SPACE_DIM>
00340 void AbstractCardiacTissue<ELEMENT_DIM,SPACE_DIM>::CalculateHaloNodesFromNodeExchange()
00341 {
00342     std::set<unsigned> halos_as_set;
00343     for (unsigned proc=0; proc<PetscTools::GetNumProcs(); proc++)
00344     {
00345         halos_as_set.insert(mNodesToReceivePerProcess[proc].begin(), mNodesToReceivePerProcess[proc].end());
00346     }
00347     mHaloNodes = std::vector<unsigned>(halos_as_set.begin(), halos_as_set.end());
00348     //PRINT_VECTOR(mHaloNodes);
00349 }
00350 
00351 
00352 template <unsigned ELEMENT_DIM,unsigned SPACE_DIM>
00353 void AbstractCardiacTissue<ELEMENT_DIM,SPACE_DIM>::SetUpHaloCells(AbstractCardiacCellFactory<ELEMENT_DIM,SPACE_DIM>* pCellFactory)
00354 {
00355     if (mExchangeHalos)
00356     {
00357         mpMesh->CalculateNodeExchange(mNodesToSendPerProcess, mNodesToReceivePerProcess);
00358         //Note that the following call will not work for a TetrahedralMesh which has no concept of halo nodes.
00359         //mpMesh->GetHaloNodeIndices( mHaloNodes );
00360         CalculateHaloNodesFromNodeExchange();
00361         unsigned num_halo_nodes = mHaloNodes.size();
00362         mHaloCellsDistributed.resize( num_halo_nodes );
00363 
00364         try
00365         {
00366             for (unsigned local_index = 0; local_index < num_halo_nodes; local_index++)
00367             {
00368                 unsigned global_index = mHaloNodes[local_index];
00369                 mHaloCellsDistributed[local_index] = pCellFactory->CreateCardiacCellForNode(global_index);
00370                 mHaloCellsDistributed[local_index]->SetUsedInTissueSimulation();
00371                 mHaloGlobalToLocalIndexMap[global_index] = local_index;
00372             }
00373 
00374             // No need to call FinaliseCellCreation() as halo node cardiac cells will
00375             // never be stimulated (their values are communicated from the process that
00376             // owns them.
00377         }
00378         catch (const Exception& e)
00379         {
00380             // Errors thrown creating cells will often be process-specific
00381             PetscTools::ReplicateException(true);
00382 
00383             // Delete cells
00384             // Should really do this for other processes too, but this is all we need
00385             // to get memory testing to pass, and leaking when we're about to die isn't
00386             // that bad!
00387             for (std::vector<AbstractCardiacCell*>::iterator cell_iterator = mHaloCellsDistributed.begin();
00388                  cell_iterator != mHaloCellsDistributed.end();
00389                  ++cell_iterator)
00390             {
00391                 delete (*cell_iterator);
00392             }
00393 
00394             throw e;
00395         }
00396         PetscTools::ReplicateException(false);
00397     }
00398 }
00399 
00400 
00401 template <unsigned ELEMENT_DIM,unsigned SPACE_DIM>
00402 void AbstractCardiacTissue<ELEMENT_DIM,SPACE_DIM>::SolveCellSystems(Vec existingSolution, double time, double nextTime, bool updateVoltage)
00403 {
00404     HeartEventHandler::BeginEvent(HeartEventHandler::SOLVE_ODES);
00405 
00406     DistributedVector dist_solution = mpDistributedVectorFactory->CreateDistributedVector(existingSolution);
00407     DistributedVector::Stripe voltage(dist_solution, 0);
00408 
00409     try
00410     {
00411         for (DistributedVector::Iterator index = dist_solution.Begin();
00412              index != dist_solution.End();
00413              ++index)
00414         {
00415             // overwrite the voltage with the input value
00416             mCellsDistributed[index.Local]->SetVoltage( voltage[index] );
00417 
00418             if(!updateVoltage)
00419             {
00420                 // solve
00421                 // Note: Voltage is not be updated. The voltage is updated in the PDE solve.
00422                 mCellsDistributed[index.Local]->ComputeExceptVoltage(time, nextTime);
00423             }
00424             else
00425             {
00426                 // solve, including updating the voltage (for the operator-splitting implementation of the monodomain solver)
00427                 mCellsDistributed[index.Local]->SolveAndUpdateState(time, nextTime);
00428                 voltage[index] = mCellsDistributed[index.Local]->GetVoltage();
00429             }
00430 
00431             // update the Iionic and stimulus caches
00432             UpdateCaches(index.Global, index.Local, nextTime);
00433         }
00434 
00435         if(updateVoltage)
00436         {
00437             dist_solution.Restore();
00438         }
00439     }
00440     catch (Exception &e)
00441     {
00442         PetscTools::ReplicateException(true);
00443         throw e;
00444     }
00445 
00446     PetscTools::ReplicateException(false);
00447     HeartEventHandler::EndEvent(HeartEventHandler::SOLVE_ODES);
00448 
00449     // Communicate new state variable values to halo nodes
00450     if (mExchangeHalos)
00451     {
00452         for ( unsigned rank_offset = 1; rank_offset < PetscTools::GetNumProcs(); rank_offset++ )
00453         {
00454             unsigned send_to      = (PetscTools::GetMyRank() + rank_offset) % (PetscTools::GetNumProcs());
00455             unsigned receive_from = (PetscTools::GetMyRank() + PetscTools::GetNumProcs()- rank_offset ) % (PetscTools::GetNumProcs());
00456 
00457             unsigned number_of_cells_to_send    = mNodesToSendPerProcess[send_to].size();
00458             unsigned number_of_cells_to_receive = mNodesToReceivePerProcess[receive_from].size();
00459 
00460             // Pack
00461             if ( number_of_cells_to_send > 0 )
00462             {
00463                 unsigned send_size = 0;
00464                 for (unsigned i=0; i<number_of_cells_to_send; i++)
00465                 {
00466                     unsigned global_cell_index = mNodesToSendPerProcess[send_to][i];
00467                     send_size += mCellsDistributed[global_cell_index - mpDistributedVectorFactory->GetLow()]->GetNumberOfStateVariables();
00468                 }
00469 
00470                 double send_data[send_size];
00471 
00472                 unsigned send_index = 0;
00473                 for (unsigned cell = 0; cell < number_of_cells_to_send; cell++)
00474                 {
00475                     unsigned global_cell_index = mNodesToSendPerProcess[send_to][cell];
00476                     AbstractCardiacCell* p_cell = mCellsDistributed[global_cell_index - mpDistributedVectorFactory->GetLow()];
00477                     std::vector<double>& cell_data = p_cell->rGetStateVariables();
00478                     const unsigned num_state_vars = p_cell->GetNumberOfStateVariables();
00479                     for (unsigned state_variable = 0; state_variable < num_state_vars; state_variable++)
00480                     {
00481                         send_data[send_index++] = cell_data[state_variable];
00482                     }
00483                 }
00484 
00485                 // Send
00486                 int ret;
00487                 ret = MPI_Send( send_data,
00488                                 send_size,
00489                                 MPI_DOUBLE,
00490                                 send_to,
00491                                 0,
00492                                 PETSC_COMM_WORLD );
00493                 assert ( ret == MPI_SUCCESS );
00494             }
00495 
00496             if ( number_of_cells_to_receive > 0 )
00497             {
00498                 // Receive
00499                 unsigned receive_size = 0;
00500                 for (unsigned i=0; i<number_of_cells_to_receive; i++)
00501                 {
00502                     unsigned halo_cell_index = mHaloGlobalToLocalIndexMap[mNodesToReceivePerProcess[receive_from][i]];
00503                     receive_size += mHaloCellsDistributed[halo_cell_index]->GetNumberOfStateVariables();
00504                 }
00505 
00506                 double receive_data[receive_size];
00507                 MPI_Status status;
00508 
00509                 int ret;
00510                 ret = MPI_Recv( receive_data,
00511                                 receive_size,
00512                                 MPI_DOUBLE,
00513                                 receive_from,
00514                                 0,
00515                                 PETSC_COMM_WORLD,
00516                                 &status );
00517                 assert ( ret == MPI_SUCCESS);
00518 
00519                 // Unpack
00520                 unsigned receive_index = 0;
00521                 for ( unsigned cell = 0; cell < number_of_cells_to_receive; cell++ )
00522                 {
00523                     AbstractCardiacCell* p_cell = mHaloCellsDistributed[mHaloGlobalToLocalIndexMap[mNodesToReceivePerProcess[receive_from][cell]]];
00524                     std::vector<double> cell_data;
00525                     cell_data.resize(p_cell->GetNumberOfStateVariables());
00526                     const unsigned number_of_state_variables = p_cell->GetNumberOfStateVariables();
00527                     for (unsigned state_variable = 0; state_variable < number_of_state_variables; state_variable++)
00528                     {
00529                         cell_data[state_variable] = receive_data[receive_index++];
00530                     }
00531                     p_cell->SetStateVariables(cell_data);
00532                 }
00533             }
00534         }
00535     }
00536 
00537     HeartEventHandler::BeginEvent(HeartEventHandler::COMMUNICATION);
00538     if ( mDoCacheReplication )
00539     {
00540         ReplicateCaches();
00541     }
00542     HeartEventHandler::EndEvent(HeartEventHandler::COMMUNICATION);
00543 }
00544 
00545 template <unsigned ELEMENT_DIM,unsigned SPACE_DIM>
00546 ReplicatableVector& AbstractCardiacTissue<ELEMENT_DIM,SPACE_DIM>::rGetIionicCacheReplicated()
00547 {
00548     return mIionicCacheReplicated;
00549 }
00550 
00551 template <unsigned ELEMENT_DIM,unsigned SPACE_DIM>
00552 ReplicatableVector& AbstractCardiacTissue<ELEMENT_DIM,SPACE_DIM>::rGetIntracellularStimulusCacheReplicated()
00553 {
00554     return mIntracellularStimulusCacheReplicated;
00555 }
00556 
00557 template <unsigned ELEMENT_DIM,unsigned SPACE_DIM>
00558 void AbstractCardiacTissue<ELEMENT_DIM,SPACE_DIM>::UpdateCaches(unsigned globalIndex, unsigned localIndex, double nextTime)
00559 {
00560     mIionicCacheReplicated[globalIndex] = mCellsDistributed[localIndex]->GetIIonic();
00561     mIntracellularStimulusCacheReplicated[globalIndex] = mCellsDistributed[localIndex]->GetIntracellularStimulus(nextTime);
00562 }
00563 
00564 template <unsigned ELEMENT_DIM,unsigned SPACE_DIM>
00565 void AbstractCardiacTissue<ELEMENT_DIM,SPACE_DIM>::ReplicateCaches()
00566 {
00567     mIionicCacheReplicated.Replicate(mpDistributedVectorFactory->GetLow(), mpDistributedVectorFactory->GetHigh());
00568     mIntracellularStimulusCacheReplicated.Replicate(mpDistributedVectorFactory->GetLow(), mpDistributedVectorFactory->GetHigh());
00569 }
00570 
00571 template <unsigned ELEMENT_DIM,unsigned SPACE_DIM>
00572 const std::vector<AbstractCardiacCell*>& AbstractCardiacTissue<ELEMENT_DIM,SPACE_DIM>::rGetCellsDistributed() const
00573 {
00574     return mCellsDistributed;
00575 }
00576 
00577 template <unsigned ELEMENT_DIM,unsigned SPACE_DIM>
00578 const AbstractTetrahedralMesh<ELEMENT_DIM,SPACE_DIM>* AbstractCardiacTissue<ELEMENT_DIM,SPACE_DIM>::pGetMesh() const
00579 {
00580     return mpMesh;
00581 }
00582 
00583 template <unsigned ELEMENT_DIM,unsigned SPACE_DIM>
00584 void AbstractCardiacTissue<ELEMENT_DIM,SPACE_DIM>::SetConductivityModifier(AbstractConductivityModifier<ELEMENT_DIM,SPACE_DIM>* pModifier)
00585 {
00586     assert(pModifier!=NULL);
00587     assert(mpConductivityModifier==NULL); // shouldn't be called twice for example, or with two different modifiers (remove this assert
00588                                           // if for whatever reason want to be able to overwrite modifiers
00589     mpConductivityModifier = pModifier;
00590 }
00591 
00592 
00594 // Explicit instantiation
00596 
00597 template class AbstractCardiacTissue<1,1>;
00598 template class AbstractCardiacTissue<1,2>;
00599 template class AbstractCardiacTissue<1,3>;
00600 template class AbstractCardiacTissue<2,2>;
00601 template class AbstractCardiacTissue<3,3>;

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