Chaste  Release::2018.1
AbstractCardiacTissue.cpp
1 /*
2 
3 Copyright (c) 2005-2018, University of Oxford.
4 All rights reserved.
5 
6 University of Oxford means the Chancellor, Masters and Scholars of the
7 University of Oxford, having an administrative office at Wellington
8 Square, Oxford OX1 2JD, UK.
9 
10 This file is part of Chaste.
11 
12 Redistribution and use in source and binary forms, with or without
13 modification, are permitted provided that the following conditions are met:
14  * Redistributions of source code must retain the above copyright notice,
15  this list of conditions and the following disclaimer.
16  * Redistributions in binary form must reproduce the above copyright notice,
17  this list of conditions and the following disclaimer in the documentation
18  and/or other materials provided with the distribution.
19  * Neither the name of the University of Oxford nor the names of its
20  contributors may be used to endorse or promote products derived from this
21  software without specific prior written permission.
22 
23 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
24 AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
25 IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
26 ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
27 LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
28 CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE
29 GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
30 HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
31 LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
32 OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
33 
34 */
35 
36 #include "AbstractCardiacTissue.hpp"
37 
38 #include <boost/scoped_array.hpp>
39 
40 #include "DistributedVector.hpp"
41 #include "AxisymmetricConductivityTensors.hpp"
42 #include "OrthotropicConductivityTensors.hpp"
43 #include "Exception.hpp"
44 #include "ChastePoint.hpp"
45 #include "AbstractChasteRegion.hpp"
46 #include "HeartEventHandler.hpp"
47 #include "PetscTools.hpp"
48 #include "PetscVecTools.hpp"
49 #include "AbstractCvodeCell.hpp"
50 #include "Warnings.hpp"
51 
52 template <unsigned ELEMENT_DIM,unsigned SPACE_DIM>
55  bool exchangeHalos)
56  : mpMesh(pCellFactory->GetMesh()),
57  mpDistributedVectorFactory(mpMesh->GetDistributedVectorFactory()),
58  mpConductivityModifier(NULL),
59  mHasPurkinje(false),
60  mDoCacheReplication(true),
61  mMeshUnarchived(false),
62  mExchangeHalos(exchangeHalos)
63 {
64  //This constructor is called from the Initialise() method of the CardiacProblem class
65  assert(pCellFactory != NULL);
66  assert(pCellFactory->GetMesh() != NULL);
67 
69  {
70  //Remove the request for a halo exchange
71  mExchangeHalos = false;
72  }
73 
74  unsigned num_local_nodes = mpDistributedVectorFactory->GetLocalOwnership();
75  bool process_has_no_nodes = (num_local_nodes == 0u);
76  if (PetscTools::ReplicateBool(process_has_no_nodes))
77  {
78  /* If there are no nodes on a process then there is a potential for an error in preconditioning.
79  * This is dangerous because the process without nodes may segfault and not propagate the error
80  * to the other processes. Therefore, to avoid deadlock, we share this potential for error between
81  * processes and throw an exception.
82  */
83 // LCOV_EXCL_START
84  // This problem normally occurs on 3 or more processes, so we can't cover it - coverage only runs with 1 and 2 processes.
85  EXCEPTION("No cells were assigned some process in AbstractCardiacTissue constructor. Advice: Make total number of processors no greater than number of nodes in the mesh");
86 // LCOV_EXCL_STOP
87  }
88  unsigned ownership_range_low = mpDistributedVectorFactory->GetLow();
89  mCellsDistributed.resize(num_local_nodes);
90 
91  // Figure out if we're dealing with Purkinje
93  = dynamic_cast<AbstractPurkinjeCellFactory<ELEMENT_DIM,SPACE_DIM>*>(pCellFactory);
94  if (p_purkinje_cell_factory)
95  {
96  mHasPurkinje = true;
97  mPurkinjeCellsDistributed.resize(num_local_nodes);
98  }
99 
101  // Set up cells
103  try
104  {
105  for (unsigned local_index = 0; local_index < num_local_nodes; local_index++)
106  {
107  unsigned global_index = ownership_range_low + local_index;
108  Node<SPACE_DIM>* p_node = mpMesh->GetNode(global_index);
109  mCellsDistributed[local_index] = pCellFactory->CreateCardiacCellForNode(p_node);
110  mCellsDistributed[local_index]->SetUsedInTissueSimulation();
111 
112  if (mHasPurkinje)
113  {
114  mPurkinjeCellsDistributed[local_index]
115  = p_purkinje_cell_factory->CreatePurkinjeCellForNode(p_node, mCellsDistributed[local_index]);
116  mPurkinjeCellsDistributed[local_index]->SetUsedInTissueSimulation();
117  }
118  }
119 
123  if (mHasPurkinje)
124  {
128  }
129  }
130  catch (const Exception& e)
131  {
132  // This catch statement is quite tricky to cover, but it is actually done
133  // in TestCardiacSimulation::TestMono1dSodiumBlockBySettingNamedParameter()
134 
135  // Errors thrown creating cells will often be process-specific
137 
138  // Delete cells
139  // Should really do this for other processes too, but this is all we need
140  // to get memory testing to pass, and leaking when we're about to die isn't
141  // that bad!
142  for (std::vector<AbstractCardiacCellInterface*>::iterator cell_iterator = mCellsDistributed.begin();
143  cell_iterator != mCellsDistributed.end();
144  ++cell_iterator)
145  {
146  delete (*cell_iterator);
147  }
148 
149  throw e;
150  }
152 
153  // Halo nodes (if required)
154  SetUpHaloCells(pCellFactory);
155 
156  HeartEventHandler::BeginEvent(HeartEventHandler::COMMUNICATION);
159 
160  if (mHasPurkinje)
161  {
164  }
165  HeartEventHandler::EndEvent(HeartEventHandler::COMMUNICATION);
166 
167  if (HeartConfig::Instance()->IsMeshProvided() && HeartConfig::Instance()->GetLoadMesh())
168  {
170  }
171  else
172  {
173  // As of r10671 fibre orientation can only be defined when loading a mesh from disc.
175  }
177 }
178 
179 // Constructor used for archiving
180 template <unsigned ELEMENT_DIM,unsigned SPACE_DIM>
182  : mpMesh(pMesh),
183  mpDistributedVectorFactory(mpMesh->GetDistributedVectorFactory()),
184  mHasPurkinje(false),
185  mDoCacheReplication(true),
186  mMeshUnarchived(true),
187  mExchangeHalos(false)
188 {
191 
194 }
195 
196 template <unsigned ELEMENT_DIM,unsigned SPACE_DIM>
198 {
199  // Delete cells
200  for (std::vector<AbstractCardiacCellInterface*>::iterator iter = mCellsDistributed.begin();
201  iter != mCellsDistributed.end();
202  ++iter)
203  {
204  delete (*iter);
205  }
206 
207  // Delete cells for halo nodes
208  for (std::vector<AbstractCardiacCellInterface*>::iterator iter = mHaloCellsDistributed.begin();
209  iter != mHaloCellsDistributed.end();
210  ++iter)
211  {
212  delete (*iter);
213  }
214 
215  delete mpIntracellularConductivityTensors;
216 
217  // Delete Purkinje cells
218  for (std::vector<AbstractCardiacCellInterface*>::iterator iter = mPurkinjeCellsDistributed.begin();
219  iter != mPurkinjeCellsDistributed.end();
220  ++iter)
221  {
222  delete (*iter);
223  }
224 
225  // If the mesh was unarchived we need to free it explicitly.
226  if (mMeshUnarchived)
227  {
228  delete mpMesh;
229  }
230 }
231 
232 
233 template <unsigned ELEMENT_DIM,unsigned SPACE_DIM>
235 {
236  return mHasPurkinje;
237 }
238 
239 template <unsigned ELEMENT_DIM,unsigned SPACE_DIM>
241 {
242  HeartEventHandler::BeginEvent(HeartEventHandler::READ_MESH);
243  mpConfig = HeartConfig::Instance();
244 
245  if (mpConfig->IsMeshProvided() && mpConfig->GetLoadMesh())
246  {
247  assert(mFibreFilePathNoExtension != "");
248 
249  switch (mpConfig->GetConductivityMedia())
250  {
251  case cp::media_type::Orthotropic:
252  {
253  mpIntracellularConductivityTensors = new OrthotropicConductivityTensors<ELEMENT_DIM,SPACE_DIM>;
254  FileFinder ortho_file(mFibreFilePathNoExtension + ".ortho", RelativeTo::AbsoluteOrCwd);
255  assert(ortho_file.Exists());
256  mpIntracellularConductivityTensors->SetFibreOrientationFile(ortho_file);
257  break;
258  }
259 
260  case cp::media_type::Axisymmetric:
261  {
262  mpIntracellularConductivityTensors = new AxisymmetricConductivityTensors<ELEMENT_DIM,SPACE_DIM>;
263  FileFinder axi_file(mFibreFilePathNoExtension + ".axi", RelativeTo::AbsoluteOrCwd);
264  assert(axi_file.Exists());
265  mpIntracellularConductivityTensors->SetFibreOrientationFile(axi_file);
266  break;
267  }
268 
269  case cp::media_type::NoFibreOrientation:
271  mpIntracellularConductivityTensors = new OrthotropicConductivityTensors<ELEMENT_DIM,SPACE_DIM>;
272  break;
273 
274  default :
276  }
277  }
278  else // Slab defined in config file or SetMesh() called; no fibre orientation assumed
279  {
281  mpIntracellularConductivityTensors = new OrthotropicConductivityTensors<ELEMENT_DIM,SPACE_DIM>;
282  }
283 
284  c_vector<double, SPACE_DIM> intra_conductivities;
285  mpConfig->GetIntracellularConductivities(intra_conductivities);
286 
287  // this definition must be here (and not inside the if statement) because SetNonConstantConductivities() will keep
288  // a pointer to it and we don't want it to go out of scope before Init() is called
289  unsigned num_local_elements = mpMesh->GetNumLocalElements();
290  std::vector<c_vector<double, SPACE_DIM> > hetero_intra_conductivities;
291 
292  if (mpConfig->GetConductivityHeterogeneitiesProvided())
293  {
294  try
295  {
296  assert(hetero_intra_conductivities.size()==0);
297  hetero_intra_conductivities.resize(num_local_elements, intra_conductivities);
298  }
299  // LCOV_EXCL_START
300  catch(std::bad_alloc &r_bad_alloc)
301  {
302  std::cout << "Failed to allocate std::vector of size " << num_local_elements << std::endl;
304  throw r_bad_alloc;
305  }
306  // LCOV_EXCL_STOP
307 
309 
310  std::vector<boost::shared_ptr<AbstractChasteRegion<SPACE_DIM> > > conductivities_heterogeneity_areas;
311  std::vector< c_vector<double,3> > intra_h_conductivities;
312  std::vector< c_vector<double,3> > extra_h_conductivities;
313  HeartConfig::Instance()->GetConductivityHeterogeneities(conductivities_heterogeneity_areas,
314  intra_h_conductivities,
315  extra_h_conductivities);
316 
317  unsigned local_element_index = 0;
318 
319  for (typename AbstractTetrahedralMesh<ELEMENT_DIM,SPACE_DIM>::ElementIterator it = mpMesh->GetElementIteratorBegin();
320  it != mpMesh->GetElementIteratorEnd();
321  ++it)
322  {
323 // unsigned element_index = it->GetIndex();
324  // if element centroid is contained in the region
325  ChastePoint<SPACE_DIM> element_centroid(it->CalculateCentroid());
326  for (unsigned region_index=0; region_index< conductivities_heterogeneity_areas.size(); region_index++)
327  {
328  if (conductivities_heterogeneity_areas[region_index]->DoesContain(element_centroid))
329  {
330  //We don't use ublas vector assignment here, because we might be getting a subvector of a 3-vector
331  for (unsigned i=0; i<SPACE_DIM; i++)
332  {
333  hetero_intra_conductivities[local_element_index][i] = intra_h_conductivities[region_index][i];
334  }
335  }
336  }
337  local_element_index++;
338  }
339 
340  mpIntracellularConductivityTensors->SetNonConstantConductivities(&hetero_intra_conductivities);
341  }
342  else
343  {
344  mpIntracellularConductivityTensors->SetConstantConductivities(intra_conductivities);
345  }
346 
347  mpIntracellularConductivityTensors->Init(this->mpMesh);
348  HeartEventHandler::EndEvent(HeartEventHandler::READ_MESH);
349 }
350 
351 
352 template <unsigned ELEMENT_DIM,unsigned SPACE_DIM>
354 {
355  mDoCacheReplication = doCacheReplication;
356 }
357 
358 template <unsigned ELEMENT_DIM,unsigned SPACE_DIM>
360 {
361  return mDoCacheReplication;
362 }
363 
364 template <unsigned ELEMENT_DIM,unsigned SPACE_DIM>
365 const c_matrix<double, SPACE_DIM, SPACE_DIM>& AbstractCardiacTissue<ELEMENT_DIM,SPACE_DIM>::rGetIntracellularConductivityTensor(unsigned elementIndex)
366 {
367  assert( mpIntracellularConductivityTensors);
368  if (mpConductivityModifier==NULL)
369  {
370  return (*mpIntracellularConductivityTensors)[elementIndex];
371  }
372  else
373  {
374  return mpConductivityModifier->rGetModifiedConductivityTensor(elementIndex, (*mpIntracellularConductivityTensors)[elementIndex], 0u);
375  }
376 }
377 
378 template <unsigned ELEMENT_DIM,unsigned SPACE_DIM>
379 const c_matrix<double, SPACE_DIM, SPACE_DIM>& AbstractCardiacTissue<ELEMENT_DIM,SPACE_DIM>::rGetExtracellularConductivityTensor(unsigned elementIndex)
380 {
381  EXCEPTION("Monodomain tissues do not have extracellular conductivity tensors.");
382 }
383 
384 template <unsigned ELEMENT_DIM,unsigned SPACE_DIM>
386 {
387  assert(mpDistributedVectorFactory->GetLow() <= globalIndex &&
388  globalIndex < mpDistributedVectorFactory->GetHigh());
389  return mCellsDistributed[globalIndex - mpDistributedVectorFactory->GetLow()];
390 }
391 
392 template <unsigned ELEMENT_DIM,unsigned SPACE_DIM>
394 {
395  assert(mpDistributedVectorFactory->GetLow() <= globalIndex &&
396  globalIndex < mpDistributedVectorFactory->GetHigh());
397  EXCEPT_IF_NOT(mHasPurkinje);
398  return mPurkinjeCellsDistributed[globalIndex - mpDistributedVectorFactory->GetLow()];
399 }
400 
401 template <unsigned ELEMENT_DIM,unsigned SPACE_DIM>
403 {
404  std::map<unsigned, unsigned>::const_iterator node_position;
405  // First search the halo
406  if ((node_position=mHaloGlobalToLocalIndexMap.find(globalIndex)) != mHaloGlobalToLocalIndexMap.end())
407  {
408  // Found a halo node
409  return mHaloCellsDistributed[node_position->second];
410  }
411  // Then search the owned node
412  if (mpDistributedVectorFactory->IsGlobalIndexLocal(globalIndex))
413  {
414  // Found an owned node
415  return mCellsDistributed[globalIndex - mpDistributedVectorFactory->GetLow()];
416  }
417  // Not here
418  EXCEPTION("Requested node/halo " << globalIndex << " does not belong to processor " << PetscTools::GetMyRank());
419 }
420 
421 
422 template <unsigned ELEMENT_DIM,unsigned SPACE_DIM>
424 {
425  std::set<unsigned> halos_as_set;
426  for (unsigned proc=0; proc<PetscTools::GetNumProcs(); proc++)
427  {
428  halos_as_set.insert(mNodesToReceivePerProcess[proc].begin(), mNodesToReceivePerProcess[proc].end());
429  }
430  mHaloNodes = std::vector<unsigned>(halos_as_set.begin(), halos_as_set.end());
431  //PRINT_VECTOR(mHaloNodes);
432 }
433 
434 
435 template <unsigned ELEMENT_DIM,unsigned SPACE_DIM>
437 {
438  if (mExchangeHalos)
439  {
440  mpMesh->CalculateNodeExchange(mNodesToSendPerProcess, mNodesToReceivePerProcess);
441  // Note that the following call will not work for a TetrahedralMesh which has
442  // no concept of halo nodes.
443  //mpMesh->GetHaloNodeIndices( mHaloNodes );
444  CalculateHaloNodesFromNodeExchange();
445  unsigned num_halo_nodes = mHaloNodes.size();
446  mHaloCellsDistributed.resize( num_halo_nodes );
447  for (unsigned local_index = 0; local_index < num_halo_nodes; local_index++)
448  {
449  unsigned global_index = mHaloNodes[local_index];
450  // These are all halo nodes, so we use the "GetNodeOrHaloNode" variety of GetNode
451  Node<SPACE_DIM>* p_node = mpMesh->GetNodeOrHaloNode(global_index);
452  mHaloCellsDistributed[local_index] = pCellFactory->CreateCardiacCellForNode(p_node);
453  mHaloCellsDistributed[local_index]->SetUsedInTissueSimulation();
454  mHaloGlobalToLocalIndexMap[global_index] = local_index;
455  }
456  // No need to call FinaliseCellCreation() as halo node cardiac cells will
457  // never be stimulated (their values are communicated from the process that
458  // owns them).
459  }
460 }
461 
462 
463 template <unsigned ELEMENT_DIM,unsigned SPACE_DIM>
464 void AbstractCardiacTissue<ELEMENT_DIM,SPACE_DIM>::SolveCellSystems(Vec existingSolution, double time, double nextTime, bool updateVoltage)
465 {
466  if (mHasPurkinje)
467  {
468  // can't do Purkinje and operator splitting
469  assert(!updateVoltage);
470  // The code below assumes Purkinje is are monodomain, so the vector has two stripes.
471  // The assert will fail the first time bidomain purkinje is coded - need to decide what
472  // ordering the three stripes (V, V_purk, phi_e) are in
473  assert(PetscVecTools::GetSize(existingSolution)==2*mpMesh->GetNumNodes());
474  }
475 
476  HeartEventHandler::BeginEvent(HeartEventHandler::SOLVE_ODES);
477 
478  DistributedVector dist_solution = mpDistributedVectorFactory->CreateDistributedVector(existingSolution);
479 
481  // Solve cell models (except purkinje cell models)
483  DistributedVector::Stripe voltage(dist_solution, 0);
484  try
485  {
486  double voltage_before_update;
487  for (DistributedVector::Iterator index = dist_solution.Begin();
488  index != dist_solution.End();
489  ++index)
490  {
491  voltage_before_update = voltage[index];
492  mCellsDistributed[index.Local]->SetVoltage( voltage_before_update );
493 
494  // Added a try-catch here to provide more output to screen when an error occurs.
496  try
497  {
498  if (!updateVoltage)
499  {
500  // solve ODE system at this node.
501  // Note: Voltage is not being updated. The voltage is updated in the PDE solve.
502 #ifndef CHASTE_CVODE
503  mCellsDistributed[index.Local]->ComputeExceptVoltage(time, nextTime);
504 #else
505  // If CVODE is enabled, and this is a CVODE cell
506  // there's a chance we can recover this by doing a reset so put the above call in a try...catch.
507  try
508  {
509  mCellsDistributed[index.Local]->ComputeExceptVoltage(time, nextTime);
510  }
511  catch (Exception &e)
512  {
513  // Try an 'emergency' reset if this is a CVODE cell.
514  // See #2594 for why we think this may be necessary.
515  if (dynamic_cast<AbstractCvodeCell*>(mCellsDistributed[index.Local]))
516  {
517  // Reset the CVODE cell, this leads to a call to CVodeReInit.
518  static_cast<AbstractCvodeCell*>(mCellsDistributed[index.Local])->ResetSolver();
519  mCellsDistributed[index.Local]->ComputeExceptVoltage(time, nextTime);
520  WARNING("Global node " << index.Global << " had an ODE solving problem in t = [" << time <<
521  ", " << nextTime << "] ms. This was fixed by a reset of CVODE, but may suggest PDE time"
522  " step should be reduced, or CVODE tolerances relaxed.");
523  }
524  else
525  {
526  throw e;
527  }
528  }
529 #endif // CHASTE_CVODE
530  }
531  else
532  {
533  // solve, including updating the voltage (for the operator-splitting implementation of the monodomain solver)
534  mCellsDistributed[index.Local]->SolveAndUpdateState(time, nextTime);
535  voltage[index] = mCellsDistributed[index.Local]->GetVoltage();
536  }
537  }
538  catch (Exception &e)
539  {
540  std::cout << std::setprecision(16);
541  std::cout << "Global node " << index.Global << " had problems with ODE solve between "
542  "t = " << time << " and " << nextTime << "ms.\n";
543 
544  std::cout << "Voltage at this node before solve was " << voltage_before_update << "mV\n"
545  "(this SHOULD NOT necessarily be the same as the one in the state variables,\n"
546  "which can be ignored and stay at the initial condition - the voltage is dictated by PDE instead of state variable.)\n";
547 
548  std::cout << "Stimulus current (NB converted to micro-Amps per cm^3) applied here is equal to:\n\t"
549  << mCellsDistributed[index.Local]->GetIntracellularStimulus(time) << " at t = " << time << "ms,\n\t"
550  << mCellsDistributed[index.Local]->GetIntracellularStimulus(nextTime) << " at t = " << nextTime << "ms.\n";
551 
552  std::cout << "Cell model: " << dynamic_cast<AbstractUntemplatedParameterisedSystem*>(mCellsDistributed[index.Local])->GetSystemName() << "\n";
553 
554  std::cout << "All state variables are now:\n";
555  std::vector<double> state_vars = mCellsDistributed[index.Local]->GetStdVecStateVariables();
556  std::vector<std::string> state_var_names = mCellsDistributed[index.Local]->rGetStateVariableNames();
557  for (unsigned i=0; i<state_vars.size(); i++)
558  {
559  std::cout << "\t" << state_var_names[i] << "\t:\t" << state_vars[i] << "\n";
560  }
561  std::cout << std::flush;
562 
563  throw e;
564  }
565  // update the Iionic and stimulus caches
566  UpdateCaches(index.Global, index.Local, nextTime);
567  }
568 
569  if (updateVoltage)
570  {
571  dist_solution.Restore();
572  }
573  }
574  catch (Exception &e)
575  {
577  throw e;
578  }
579 
581  // Solve purkinje cell models
583  if (mHasPurkinje)
584  {
585  DistributedVector::Stripe purkinje_voltage(dist_solution, 1);
586  try
587  {
588  for (DistributedVector::Iterator index = dist_solution.Begin();
589  index != dist_solution.End();
590  ++index)
591  {
592  // overwrite the voltage with the input value
593  mPurkinjeCellsDistributed[index.Local]->SetVoltage( purkinje_voltage[index] );
594 
595  // solve
596  // Note: Voltage is not being updated. The voltage is updated in the PDE solve.
597  mPurkinjeCellsDistributed[index.Local]->ComputeExceptVoltage(time, nextTime);
598 
599  // update the Iionic and stimulus caches
600  UpdatePurkinjeCaches(index.Global, index.Local, nextTime);
601  }
602  }
603  // LCOV_EXCL_START
604  catch (Exception& e)
605  {
609 
612  throw e;
613  }
614  // LCOV_EXCL_STOP
615  }
616 
618  HeartEventHandler::EndEvent(HeartEventHandler::SOLVE_ODES);
619 
620  // Communicate new state variable values to halo nodes
621  if (mExchangeHalos)
622  {
623  assert(!mHasPurkinje);
624 
625  for ( unsigned rank_offset = 1; rank_offset < PetscTools::GetNumProcs(); rank_offset++ )
626  {
627  unsigned send_to = (PetscTools::GetMyRank() + rank_offset) % (PetscTools::GetNumProcs());
628  unsigned receive_from = (PetscTools::GetMyRank() + PetscTools::GetNumProcs()- rank_offset ) % (PetscTools::GetNumProcs());
629 
630  unsigned number_of_cells_to_send = mNodesToSendPerProcess[send_to].size();
631  unsigned number_of_cells_to_receive = mNodesToReceivePerProcess[receive_from].size();
632 
633  // Pack send buffer
634  unsigned send_size = 0;
635  for (unsigned i=0; i<number_of_cells_to_send; i++)
636  {
637  unsigned global_cell_index = mNodesToSendPerProcess[send_to][i];
638  send_size += mCellsDistributed[global_cell_index - mpDistributedVectorFactory->GetLow()]->GetNumberOfStateVariables();
639  }
640 
641  boost::scoped_array<double> send_data(new double[send_size]);
642 
643  unsigned send_index = 0;
644  for (unsigned cell = 0; cell < number_of_cells_to_send; cell++)
645  {
646  unsigned global_cell_index = mNodesToSendPerProcess[send_to][cell];
647  AbstractCardiacCellInterface* p_cell = mCellsDistributed[global_cell_index - mpDistributedVectorFactory->GetLow()];
648  std::vector<double> cell_data = p_cell->GetStdVecStateVariables();
649  const unsigned num_state_vars = p_cell->GetNumberOfStateVariables();
650  for (unsigned state_variable = 0; state_variable < num_state_vars; state_variable++)
651  {
652  send_data[send_index++] = cell_data[state_variable];
653  }
654  }
655  // Receive buffer
656  unsigned receive_size = 0;
657  for (unsigned i=0; i<number_of_cells_to_receive; i++)
658  {
659  unsigned halo_cell_index = mHaloGlobalToLocalIndexMap[mNodesToReceivePerProcess[receive_from][i]];
660  receive_size += mHaloCellsDistributed[halo_cell_index]->GetNumberOfStateVariables();
661  }
662 
663  boost::scoped_array<double> receive_data(new double[receive_size]);
664 
665  // Send and receive
666  int ret;
667  MPI_Status status;
668  ret = MPI_Sendrecv(send_data.get(), send_size,
669  MPI_DOUBLE,
670  send_to, 0,
671  receive_data.get(), receive_size,
672  MPI_DOUBLE,
673  receive_from, 0,
674  PETSC_COMM_WORLD, &status);
675  UNUSED_OPT(ret);
676  assert ( ret == MPI_SUCCESS);
677 
678  // Unpack
679  unsigned receive_index = 0;
680  for ( unsigned cell = 0; cell < number_of_cells_to_receive; cell++ )
681  {
682  AbstractCardiacCellInterface* p_cell = mHaloCellsDistributed[mHaloGlobalToLocalIndexMap[mNodesToReceivePerProcess[receive_from][cell]]];
683  const unsigned number_of_state_variables = p_cell->GetNumberOfStateVariables();
684 
685  std::vector<double> cell_data(number_of_state_variables);
686  for (unsigned state_variable = 0; state_variable < number_of_state_variables; state_variable++)
687  {
688  cell_data[state_variable] = receive_data[receive_index++];
689  }
690  p_cell->SetStateVariables(cell_data);
691  }
692  }
693  }
694 
695  HeartEventHandler::BeginEvent(HeartEventHandler::COMMUNICATION);
696  if (mDoCacheReplication)
697  {
698  ReplicateCaches();
699  }
700  HeartEventHandler::EndEvent(HeartEventHandler::COMMUNICATION);
701 }
702 
703 template <unsigned ELEMENT_DIM,unsigned SPACE_DIM>
705 {
706  return mIionicCacheReplicated;
707 }
708 
709 template <unsigned ELEMENT_DIM,unsigned SPACE_DIM>
711 {
712  return mIntracellularStimulusCacheReplicated;
713 }
714 
715 template <unsigned ELEMENT_DIM,unsigned SPACE_DIM>
717 {
718  EXCEPT_IF_NOT(mHasPurkinje);
719  return mPurkinjeIionicCacheReplicated;
720 }
721 
722 template <unsigned ELEMENT_DIM,unsigned SPACE_DIM>
724 {
725  EXCEPT_IF_NOT(mHasPurkinje);
726  return mPurkinjeIntracellularStimulusCacheReplicated;
727 }
728 
729 template <unsigned ELEMENT_DIM,unsigned SPACE_DIM>
730 void AbstractCardiacTissue<ELEMENT_DIM,SPACE_DIM>::UpdateCaches(unsigned globalIndex, unsigned localIndex, double nextTime)
731 {
732  mIionicCacheReplicated[globalIndex] = mCellsDistributed[localIndex]->GetIIonic();
733  mIntracellularStimulusCacheReplicated[globalIndex] = mCellsDistributed[localIndex]->GetIntracellularStimulus(nextTime);
734 }
735 
736 template <unsigned ELEMENT_DIM,unsigned SPACE_DIM>
737 void AbstractCardiacTissue<ELEMENT_DIM,SPACE_DIM>::UpdatePurkinjeCaches(unsigned globalIndex, unsigned localIndex, double nextTime)
738 {
739  assert(mHasPurkinje);
740  mPurkinjeIionicCacheReplicated[globalIndex] = mPurkinjeCellsDistributed[localIndex]->GetIIonic();
741  mPurkinjeIntracellularStimulusCacheReplicated[globalIndex] = mPurkinjeCellsDistributed[localIndex]->GetIntracellularStimulus(nextTime);
742 }
743 
744 template <unsigned ELEMENT_DIM,unsigned SPACE_DIM>
746 {
747  // ReplicateCaches only needed for SVI (and non-matrix based assembly which is no longer in code)
748  // which is not implemented with Purkinje. See commented code below if introducing this.
749  assert(!mHasPurkinje);
750 
751  mIionicCacheReplicated.Replicate(mpDistributedVectorFactory->GetLow(), mpDistributedVectorFactory->GetHigh());
752  mIntracellularStimulusCacheReplicated.Replicate(mpDistributedVectorFactory->GetLow(), mpDistributedVectorFactory->GetHigh());
753 
754  //if (mHasPurkinje)
755  //{
756  // mPurkinjeIionicCacheReplicated.Replicate(mpDistributedVectorFactory->GetLow(), mpDistributedVectorFactory->GetHigh());
757  // mPurkinjeIntracellularStimulusCacheReplicated.Replicate(mpDistributedVectorFactory->GetLow(), mpDistributedVectorFactory->GetHigh());
758  //}
759 }
760 
761 template <unsigned ELEMENT_DIM,unsigned SPACE_DIM>
762 const std::vector<AbstractCardiacCellInterface*>& AbstractCardiacTissue<ELEMENT_DIM,SPACE_DIM>::rGetCellsDistributed() const
763 {
764  return mCellsDistributed;
765 }
766 
767 template <unsigned ELEMENT_DIM,unsigned SPACE_DIM>
768 const std::vector<AbstractCardiacCellInterface*>& AbstractCardiacTissue<ELEMENT_DIM,SPACE_DIM>::rGetPurkinjeCellsDistributed() const
769 {
770  EXCEPT_IF_NOT(mHasPurkinje);
771  return mPurkinjeCellsDistributed;
772 }
773 
774 template <unsigned ELEMENT_DIM,unsigned SPACE_DIM>
776 {
777  return mpMesh;
778 }
779 
780 template <unsigned ELEMENT_DIM,unsigned SPACE_DIM>
782 {
783  assert(pModifier!=NULL);
784  assert(mpConductivityModifier==NULL); // shouldn't be called twice for example, or with two different modifiers (remove this assert
785  // if for whatever reason want to be able to overwrite modifiers)
786  mpConductivityModifier = pModifier;
787 }
788 
789 // Explicit instantiation
790 template class AbstractCardiacTissue<1,1>;
791 template class AbstractCardiacTissue<1,2>;
792 template class AbstractCardiacTissue<1,3>;
793 template class AbstractCardiacTissue<2,2>;
794 template class AbstractCardiacTissue<3,3>;
virtual void FinaliseCellCreation(std::vector< AbstractCardiacCellInterface * > *pCellsDistributed, unsigned lo, unsigned hi)
ReplicatableVector mIntracellularStimulusCacheReplicated
virtual AbstractCardiacCellInterface * CreateCardiacCellForNode(Node< SPACE_DIM > *pNode)
static unsigned GetSize(Vec vector)
void SetConductivityModifier(AbstractConductivityModifier< ELEMENT_DIM, SPACE_DIM > *pModifier)
DistributedVectorFactory * mpDistributedVectorFactory
AbstractCardiacCellInterface * CreatePurkinjeCellForNode(Node< SPACE_DIM > *pNode, AbstractCardiacCellInterface *pCardiacCell)
Definition: Node.hpp:58
ReplicatableVector & rGetPurkinjeIionicCacheReplicated()
static bool ReplicateBool(bool flag)
Definition: PetscTools.cpp:186
AbstractCardiacTissue(AbstractCardiacCellFactory< ELEMENT_DIM, SPACE_DIM > *pCellFactory, bool exchangeHalos=false)
AbstractTetrahedralMesh< ELEMENT_DIM, SPACE_DIM > * GetMesh()
virtual void SolveCellSystems(Vec existingSolution, double time, double nextTime, bool updateVoltage=false)
static std::string GetMeshFilename()
#define EXCEPTION(message)
Definition: Exception.hpp:143
void UpdateCaches(unsigned globalIndex, unsigned localIndex, double nextTime)
ReplicatableVector mPurkinjeIionicCacheReplicated
std::vector< AbstractCardiacCellInterface * > mPurkinjeCellsDistributed
virtual std::vector< double > GetStdVecStateVariables()=0
const std::vector< std::string > & rGetStateVariableNames() const
#define NEVER_REACHED
Definition: Exception.hpp:206
ReplicatableVector mPurkinjeIntracellularStimulusCacheReplicated
std::string GetMeshName() const
void GetConductivityHeterogeneities(std::vector< boost::shared_ptr< AbstractChasteRegion< DIM > > > &conductivitiesHeterogeneityAreas, std::vector< c_vector< double, 3 > > &intraConductivities, std::vector< c_vector< double, 3 > > &extraConductivities) const
virtual void SetStateVariables(const std::vector< double > &rVariables)=0
virtual unsigned GetNumberOfStateVariables() const =0
std::vector< AbstractCardiacCellInterface * > mCellsDistributed
static bool IsSequential()
Definition: PetscTools.cpp:91
AbstractCardiacCellInterface * GetCardiacCellOrHaloCell(unsigned globalIndex)
ReplicatableVector & rGetIntracellularStimulusCacheReplicated()
AbstractCardiacCellInterface * GetPurkinjeCell(unsigned globalIndex)
void SetUpHaloCells(AbstractCardiacCellFactory< ELEMENT_DIM, SPACE_DIM > *pCellFactory)
virtual void FinalisePurkinjeCellCreation(std::vector< AbstractCardiacCellInterface * > *pPurkinjeCellsDistributed, unsigned lo, unsigned hi)
static void ReplicateException(bool flag)
Definition: PetscTools.cpp:198
ReplicatableVector mIionicCacheReplicated
AbstractTetrahedralMesh< ELEMENT_DIM, SPACE_DIM > * mpMesh
void ComputeExceptVoltage(double tStart, double tEnd)
#define EXCEPT_IF_NOT(test)
Definition: Exception.hpp:158
const std::vector< AbstractCardiacCellInterface * > & rGetCellsDistributed() const
virtual const c_matrix< double, SPACE_DIM, SPACE_DIM > & rGetExtracellularConductivityTensor(unsigned elementIndex)
bool Exists() const
Definition: FileFinder.cpp:183
#define UNUSED_OPT(var)
Definition: Exception.hpp:216
void Resize(unsigned size)
const std::vector< AbstractCardiacCellInterface * > & rGetPurkinjeCellsDistributed() const
void SetCacheReplication(bool doCacheReplication)
ReplicatableVector & rGetPurkinjeIntracellularStimulusCacheReplicated()
static unsigned GetNumProcs()
Definition: PetscTools.cpp:108
const c_matrix< double, SPACE_DIM, SPACE_DIM > & rGetIntracellularConductivityTensor(unsigned elementIndex)
static unsigned GetMyRank()
Definition: PetscTools.cpp:114
static std::string GetArchiveDirectory()
AbstractCardiacCellInterface * GetCardiacCell(unsigned globalIndex)
void UpdatePurkinjeCaches(unsigned globalIndex, unsigned localIndex, double nextTime)
static HeartConfig * Instance()
const AbstractTetrahedralMesh< ELEMENT_DIM, SPACE_DIM > * pGetMesh() const
ReplicatableVector & rGetIionicCacheReplicated()