Chaste  Release::3.4
CardiacElectroMechanicsProblem.cpp
1 /*
2 
3 Copyright (c) 2005-2016, 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 "CardiacElectroMechanicsProblem.hpp"
37 
38 #include "OutputFileHandler.hpp"
39 #include "ReplicatableVector.hpp"
40 #include "HeartConfig.hpp"
41 #include "LogFile.hpp"
42 #include "ChastePoint.hpp"
43 #include "Element.hpp"
44 #include "BoundaryConditionsContainer.hpp"
45 #include "AbstractDynamicLinearPdeSolver.hpp"
46 #include "TimeStepper.hpp"
47 #include "TrianglesMeshWriter.hpp"
48 #include "Hdf5ToMeshalyzerConverter.hpp"
49 #include "Hdf5ToCmguiConverter.hpp"
50 #include "MeshalyzerMeshWriter.hpp"
51 #include "PetscTools.hpp"
52 #include "ImplicitCardiacMechanicsSolver.hpp"
53 #include "ExplicitCardiacMechanicsSolver.hpp"
54 #include "CmguiDeformedSolutionsWriter.hpp"
55 #include "VoltageInterpolaterOntoMechanicsMesh.hpp"
56 #include "BidomainProblem.hpp"
57 
58 
59 template<unsigned DIM, unsigned ELEC_PROB_DIM>
61 {
62  assert(mIsWatchedLocation);
63 
64  // find the nearest electrics mesh node
65  double min_dist = DBL_MAX;
66  unsigned node_index = UNSIGNED_UNSET;
67  for(unsigned i=0; i<mpElectricsMesh->GetNumNodes(); i++)
68  {
69  double dist = norm_2(mWatchedLocation - mpElectricsMesh->GetNode(i)->rGetLocation());
70  if(dist < min_dist)
71  {
72  min_dist = dist;
73  node_index = i;
74  }
75  }
76 
77  // set up watched node, if close enough
78  assert(node_index != UNSIGNED_UNSET); // should def have found something
79  c_vector<double,DIM> pos = mpElectricsMesh->GetNode(node_index)->rGetLocation();
80 
81  if(min_dist > 1e-8)
82  {
83  #define COVERAGE_IGNORE
84  std::cout << "ERROR: Could not find an electrics node very close to requested watched location - "
85  << "min distance was " << min_dist << " for node " << node_index
86  << " at location " << pos << std::flush;;
87 
89  //EXCEPTION("Could not find an electrics node very close to requested watched location");
91  #undef COVERAGE_IGNORE
92  }
93  else
94  {
95  LOG(2,"Chosen electrics node "<<node_index<<" at location " << pos << " to be watched");
96  mWatchedElectricsNodeIndex = node_index;
97  }
98 
99  // find nearest mechanics mesh
100  min_dist = DBL_MAX;
101  node_index = UNSIGNED_UNSET;
102  c_vector<double,DIM> pos_at_min;
103 
104  for(unsigned i=0; i<mpMechanicsMesh->GetNumNodes(); i++)
105  {
106  c_vector<double,DIM> position = mpMechanicsMesh->GetNode(i)->rGetLocation();
107 
108  double dist = norm_2(position-mWatchedLocation);
109 
110  if(dist < min_dist)
111  {
112  min_dist = dist;
113  node_index = i;
114  pos_at_min = position;
115  }
116  }
117 
118  // set up watched node, if close enough
119  assert(node_index != UNSIGNED_UNSET); // should def have found something
120 
121  if(min_dist > 1e-8)
122  {
123  #define COVERAGE_IGNORE
124  std::cout << "ERROR: Could not find a mechanics node very close to requested watched location - "
125  << "min distance was " << min_dist << " for node " << node_index
126  << " at location " << pos_at_min;
127 
129  //EXCEPTION("Could not find a mechanics node very close to requested watched location");
131  #undef COVERAGE_IGNORE
132  }
133  else
134  {
135  LOG(2,"Chosen mechanics node "<<node_index<<" at location " << pos << " to be watched");
136  mWatchedMechanicsNodeIndex = node_index;
137  }
138 
139  OutputFileHandler handler(mOutputDirectory,false);
140  mpWatchedLocationFile = handler.OpenOutputFile("watched.txt");
141 }
142 
143 template<unsigned DIM, unsigned ELEC_PROB_DIM>
145 {
146  assert(mIsWatchedLocation);
147 
148  std::vector<c_vector<double,DIM> >& deformed_position = mpMechanicsSolver->rGetDeformedPosition();
149 
151  ReplicatableVector voltage_replicated(voltage);
152  double V = voltage_replicated[mWatchedElectricsNodeIndex];
153 
156  // // Metadata is currently being added to CellML models and then this will be avoided by asking for Calcium.
157  // double Ca = mpElectricsProblem->GetMonodomainTissue()->GetCardiacCell(mWatchedElectricsNodeIndex)->GetIntracellularCalciumConcentration();
158 
159  *mpWatchedLocationFile << time << " ";
160  for(unsigned i=0; i<DIM; i++)
161  {
162  *mpWatchedLocationFile << deformed_position[mWatchedMechanicsNodeIndex](i) << " ";
163  }
164  *mpWatchedLocationFile << V << "\n";
165  mpWatchedLocationFile->flush();
166 }
167 
168 template<unsigned DIM, unsigned ELEC_PROB_DIM>
169 c_matrix<double,DIM,DIM>& CardiacElectroMechanicsProblem<DIM,ELEC_PROB_DIM>::rCalculateModifiedConductivityTensor(unsigned elementIndex, const c_matrix<double,DIM,DIM>& rOriginalConductivity, unsigned domainIndex)
170 {
171 
172  // first get the deformation gradient for this electrics element
173  unsigned containing_mechanics_elem = mpMeshPair->rGetCoarseElementsForFineElementCentroids()[elementIndex];
174  c_matrix<double,DIM,DIM>& r_deformation_gradient = mDeformationGradientsForEachMechanicsElement[containing_mechanics_elem];
175 
176  // compute sigma_def = F^{-1} sigma_undef F^{-T}
177  c_matrix<double,DIM,DIM> inv_F = Inverse(r_deformation_gradient);
178  mModifiedConductivityTensor = prod(inv_F, rOriginalConductivity);
179  mModifiedConductivityTensor = prod(mModifiedConductivityTensor, trans(inv_F));
180 
181  return mModifiedConductivityTensor;
182 }
183 
184 
188 template<unsigned DIM, unsigned PROBLEM_DIM>
190 {
191 public:
198  static AbstractCardiacProblem<DIM, DIM, PROBLEM_DIM>* Create(ElectricsProblemType problemType,
199  AbstractCardiacCellFactory<DIM>* pCellFactory);
200 };
201 
205 template<unsigned DIM>
207 {
208 public:
215  static AbstractCardiacProblem<DIM, DIM, 1u>* Create(ElectricsProblemType problemType,
216  AbstractCardiacCellFactory<DIM>* pCellFactory)
217  {
218  if (problemType == MONODOMAIN)
219  {
220  return new MonodomainProblem<DIM>(pCellFactory);
221  }
222  EXCEPTION("The second template parameter should be 2 when a bidomain problem is chosen");
223  }
224 };
225 
229 template<unsigned DIM>
231 {
232 public:
239  static AbstractCardiacProblem<DIM, DIM, 2u>* Create(ElectricsProblemType problemType,
240  AbstractCardiacCellFactory<DIM>* pCellFactory)
241  {
242  if (problemType == BIDOMAIN)
243  {
244  return new BidomainProblem<DIM>(pCellFactory, false);//false-> no bath
245  }
246  if (problemType == BIDOMAIN_WITH_BATH)
247  {
248  return new BidomainProblem<DIM>(pCellFactory, true);// true-> bath
249  }
250  EXCEPTION("The second template parameter should be 1 when a monodomain problem is chosen");
251  }
252 };
253 
254 
255 template<unsigned DIM, unsigned ELEC_PROB_DIM>
257  CompressibilityType compressibilityType,
258  ElectricsProblemType electricsProblemType,
259  TetrahedralMesh<DIM,DIM>* pElectricsMesh,
260  QuadraticMesh<DIM>* pMechanicsMesh,
261  AbstractCardiacCellFactory<DIM>* pCellFactory,
262  ElectroMechanicsProblemDefinition<DIM>* pProblemDefinition,
263  std::string outputDirectory)
264  : mCompressibilityType(compressibilityType),
265  mpCardiacMechSolver(NULL),
266  mpMechanicsSolver(NULL),
267  mpElectricsMesh(pElectricsMesh),
268  mpMechanicsMesh(pMechanicsMesh),
269  mpProblemDefinition(pProblemDefinition),
270  mHasBath(false),
271  mpMeshPair(NULL),
272  mNoElectricsOutput(false),
273  mIsWatchedLocation(false),
274  mWatchedElectricsNodeIndex(UNSIGNED_UNSET),
275  mWatchedMechanicsNodeIndex(UNSIGNED_UNSET),
276  mNumTimestepsToOutputDeformationGradientsAndStress(UNSIGNED_UNSET)
277 {
278  // Do some initial set up...
279  // However, NOTE, we don't use either the passed in meshes or the problem_definition.
280  // These pointers are allowed to be NULL, in case a child constructor wants to set
281  // them up (eg CardiacElectroMechProbRegularGeom).
282  // The meshes and problem_defn are used for the first time in Initialise().
283 
284 
285  // Start-up mechanics event handler..
287  MechanicsEventHandler::BeginEvent(MechanicsEventHandler::ALL);
288  // disable the electric event handler, because we use a problem class but
289  // don't call Solve, so we would have to worry about starting and ending any
290  // events in AbstractCardiacProblem::Solve() (esp. calling EndEvent(EVERYTHING))
291  // if we didn't disable it.
293 
294  assert(HeartConfig::Instance()->GetSimulationDuration()>0.0);
295  assert(HeartConfig::Instance()->GetPdeTimeStep()>0.0);
296 
297  // Create the monodomain problem.
298  // **NOTE** WE ONLY USE THIS TO: set up the cells, get an initial condition
299  // (voltage) vector, and get a solver. We won't ever call Solve on the cardiac problem class
300  assert(pCellFactory != NULL);
301  mpElectricsProblem = CreateElectricsProblem<DIM,ELEC_PROB_DIM>::Create(electricsProblemType, pCellFactory);
302 
303  if (electricsProblemType == BIDOMAIN_WITH_BATH)
304  {
305  mHasBath = true;
306  }
307  // check whether output is required
308  mWriteOutput = (outputDirectory!="");
309  if(mWriteOutput)
310  {
311  mOutputDirectory = outputDirectory;
312  // create the directory
317  }
318  else
319  {
321  }
322 
323 // mpImpactRegion=NULL;
324 }
325 
326 template<unsigned DIM, unsigned ELEC_PROB_DIM>
328 {
329  // NOTE if SetWatchedLocation but not Initialise has been called, mpWatchedLocationFile
330  // will be uninitialised and using it will cause a seg fault. Hence the mpMechanicsMesh!=NULL
331  // it is true if Initialise has been called.
332  if(mIsWatchedLocation && mpMechanicsMesh)
333  {
334  mpWatchedLocationFile->close();
335  }
336 
337  delete mpElectricsProblem;
338  delete mpCardiacMechSolver;
339  delete mpMeshPair;
340 
341  LogFile::Close();
342 }
343 
344 template<unsigned DIM, unsigned ELEC_PROB_DIM>
346 {
347  assert(mpElectricsMesh!=NULL);
348  assert(mpMechanicsMesh!=NULL);
349  assert(mpProblemDefinition!=NULL);
350  assert(mpCardiacMechSolver==NULL);
351 
354  mNumElecTimestepsPerMechTimestep = (unsigned) floor((mpProblemDefinition->GetMechanicsSolveTimestep()/HeartConfig::Instance()->GetPdeTimeStep())+0.5);
355  if(fabs(mNumElecTimestepsPerMechTimestep*HeartConfig::Instance()->GetPdeTimeStep() - mpProblemDefinition->GetMechanicsSolveTimestep()) > 1e-6)
356  {
357  EXCEPTION("Electrics PDE timestep does not divide mechanics solve timestep");
358  }
359 
360  // Create the Logfile (note we have to do this after the output dir has been
361  // created, else the log file might get cleaned away
362  std::string log_dir = mOutputDirectory; // just the TESTOUTPUT dir if mOutputDir="";
363  LogFile::Instance()->Set(2, mOutputDirectory);
364  LogFile::Instance()->WriteHeader("Electromechanics");
365  LOG(2, DIM << "d Implicit CardiacElectroMechanics Simulation:");
366  LOG(2, "End time = " << HeartConfig::Instance()->GetSimulationDuration() << ", electrics time step = " << HeartConfig::Instance()->GetPdeTimeStep() << ", mechanics timestep = " << mpProblemDefinition->GetMechanicsSolveTimestep() << "\n");
367  LOG(2, "Contraction model ode timestep " << mpProblemDefinition->GetContractionModelOdeTimestep());
368  LOG(2, "Output is written to " << mOutputDirectory << "/[deformation/electrics]");
369 
370  LOG(2, "Electrics mesh has " << mpElectricsMesh->GetNumNodes() << " nodes");
371  LOG(2, "Mechanics mesh has " << mpMechanicsMesh->GetNumNodes() << " nodes");
372 
373  LOG(2, "Initialising..");
374 
375 
376  if(mIsWatchedLocation)
377  {
378  DetermineWatchedNodes();
379  }
380 
381  // initialise electrics problem
382  mpElectricsProblem->SetMesh(mpElectricsMesh);
383  mpElectricsProblem->Initialise();
384 
385  if(mCompressibilityType==INCOMPRESSIBLE)
386  {
387  switch(mpProblemDefinition->GetSolverType())
388  {
389  case EXPLICIT:
391  *mpMechanicsMesh,*mpProblemDefinition,mDeformationOutputDirectory);
392  break;
393  case IMPLICIT:
395  *mpMechanicsMesh,*mpProblemDefinition,mDeformationOutputDirectory);
396  break;
397  default:
399  }
400  }
401  else
402  {
403  // repeat above with Compressible solver rather than incompressible -
404  // not the neatest but avoids having to template this class.
405  switch(mpProblemDefinition->GetSolverType())
406  {
407  case EXPLICIT:
409  *mpMechanicsMesh,*mpProblemDefinition,mDeformationOutputDirectory);
410  break;
411  case IMPLICIT:
413  *mpMechanicsMesh,*mpProblemDefinition,mDeformationOutputDirectory);
414  break;
415  default:
417  }
418  }
419 
420 
421  mpMechanicsSolver = dynamic_cast<AbstractNonlinearElasticitySolver<DIM>*>(mpCardiacMechSolver);
422  assert(mpMechanicsSolver);
423 
424  // set up mesh pair and determine the fine mesh elements and corresponding weights for each
425  // quadrature point in the coarse mesh
426  mpMeshPair = new FineCoarseMeshPair<DIM>(*mpElectricsMesh, *mpMechanicsMesh);
427  mpMeshPair->SetUpBoxesOnFineMesh();
428  mpMeshPair->ComputeFineElementsAndWeightsForCoarseQuadPoints(*(mpCardiacMechSolver->GetQuadratureRule()), false);
429  mpMeshPair->DeleteFineBoxCollection();
430 
431  mpCardiacMechSolver->SetFineCoarseMeshPair(mpMeshPair);
432  mpCardiacMechSolver->Initialise();
433 
434  unsigned num_quad_points = mpCardiacMechSolver->GetTotalNumQuadPoints();
435  mInterpolatedCalciumConcs.assign(num_quad_points, 0.0);
436  mInterpolatedVoltages.assign(num_quad_points, 0.0);
437 
438  if(mpProblemDefinition->ReadFibreSheetDirectionsFromFile())
439  {
440  mpCardiacMechSolver->SetVariableFibreSheetDirections(mpProblemDefinition->GetFibreSheetDirectionsFile(),
441  mpProblemDefinition->GetFibreSheetDirectionsDefinedPerQuadraturePoint());
442  }
443 
444 
445  if(mpProblemDefinition->GetDeformationAffectsConductivity() || mpProblemDefinition->GetDeformationAffectsCellModels())
446  {
447  mpMeshPair->SetUpBoxesOnCoarseMesh();
448  }
449 
450 
451  if(mpProblemDefinition->GetDeformationAffectsCellModels() || mpProblemDefinition->GetDeformationAffectsConductivity())
452  {
453  // initialise the stretches saved for each mechanics element
454  mStretchesForEachMechanicsElement.resize(mpMechanicsMesh->GetNumElements(), 1.0);
455 
456  // initialise the store of the F in each mechanics element (one constant value of F) in each
457  mDeformationGradientsForEachMechanicsElement.resize(mpMechanicsMesh->GetNumElements(),identity_matrix<double>(DIM));
458  }
459 
460 
461  if(mpProblemDefinition->GetDeformationAffectsCellModels())
462  {
463  // compute the coarse elements which contain each fine node -- for transferring stretch from
464  // mechanics solve electrics cell models
465  mpMeshPair->ComputeCoarseElementsForFineNodes(false);
466 
467  }
468 
469  if(mpProblemDefinition->GetDeformationAffectsConductivity())
470  {
471  // compute the coarse elements which contain each fine element centroid -- for transferring F from
472  // mechanics solve to electrics mesh elements
473  mpMeshPair->ComputeCoarseElementsForFineElementCentroids(false);
474 
475  // tell the abstract tissue class that the conductivities need to be modified, passing in this class
476  // (which is of type AbstractConductivityModifier)
477  mpElectricsProblem->GetTissue()->SetConductivityModifier(this);
478  }
479 
480  if(mWriteOutput)
481  {
482  TrianglesMeshWriter<DIM,DIM> mesh_writer(mOutputDirectory,"electrics_mesh",false);
483  mesh_writer.WriteFilesUsingMesh(*mpElectricsMesh);
484  }
485 }
486 
487 template<unsigned DIM, unsigned ELEC_PROB_DIM>
489 {
490  // initialise the meshes and mechanics solver
491  if(mpCardiacMechSolver==NULL)
492  {
493  Initialise();
494  }
495 
496  bool verbose_during_solve = ( mpProblemDefinition->GetVerboseDuringSolve()
497  || CommandLineArguments::Instance()->OptionExists("-mech_verbose")
498  || CommandLineArguments::Instance()->OptionExists("-mech_very_verbose"));
499 
500 
501  mpProblemDefinition->Validate();
502 
503  boost::shared_ptr<BoundaryConditionsContainer<DIM,DIM,ELEC_PROB_DIM> > p_bcc(new BoundaryConditionsContainer<DIM,DIM,ELEC_PROB_DIM>);
504  p_bcc->DefineZeroNeumannOnMeshBoundary(mpElectricsMesh, 0);
505  mpElectricsProblem->SetBoundaryConditionsContainer(p_bcc);
506 
507  // get an electrics solver from the problem. Note that we don't call
508  // Solve() on the CardiacProblem class, we do the looping here.
510  = mpElectricsProblem->CreateSolver();
511 
512  // set up initial voltage etc
513  Vec electrics_solution=NULL; //This will be set and used later
514  Vec calcium_data= mpElectricsMesh->GetDistributedVectorFactory()->CreateVec();
515  Vec initial_voltage = mpElectricsProblem->CreateInitialCondition();
516 
517  // write the initial position
518  unsigned counter = 0;
519 
520  TimeStepper stepper(0.0, HeartConfig::Instance()->GetSimulationDuration(), mpProblemDefinition->GetMechanicsSolveTimestep());
521 
522  CmguiDeformedSolutionsWriter<DIM>* p_cmgui_writer = NULL;
523 
524  std::vector<std::string> variable_names;
525 
526  if (mWriteOutput)
527  {
528  mpMechanicsSolver->SetWriteOutput();
529  mpMechanicsSolver->WriteCurrentSpatialSolution("undeformed","nodes");
530 
531  p_cmgui_writer = new CmguiDeformedSolutionsWriter<DIM>(mOutputDirectory+"/deformation/cmgui",
532  "solution",
533  *(this->mpMechanicsMesh),
534  WRITE_QUADRATIC_MESH);
535  variable_names.push_back("V");
536  if(ELEC_PROB_DIM==2)
537  {
538  variable_names.push_back("Phi_e");
539  if (mHasBath==true)
540  {
541  std::vector<std::string> regions;
542  regions.push_back("tissue");
543  regions.push_back("bath");
544  p_cmgui_writer->SetRegionNames(regions);
545  }
546  }
547  p_cmgui_writer->SetAdditionalFieldNames(variable_names);
548  p_cmgui_writer->WriteInitialMesh("undeformed");
549  }
550 
557 
558  LOG(2, "\nSolving for initial deformation");
559  #define COVERAGE_IGNORE
560  if(verbose_during_solve)
561  {
562  std::cout << "\n\n ** Solving for initial deformation\n";
563  }
564  #undef COVERAGE_IGNORE
565 
566  mpMechanicsSolver->SetWriteOutput(false);
567 
568  mpMechanicsSolver->SetCurrentTime(0.0);
569  MechanicsEventHandler::BeginEvent(MechanicsEventHandler::ALL_MECH);
570 
571  mpMechanicsSolver->SetIncludeActiveTension(false);
572  if(mNumTimestepsToOutputDeformationGradientsAndStress!=UNSIGNED_UNSET)
573  {
574  mpMechanicsSolver->SetComputeAverageStressPerElementDuringSolve(true);
575  }
576 
577  unsigned total_newton_iters = 0;
578  for(unsigned index=1; index<=mpProblemDefinition->GetNumIncrementsForInitialDeformation(); index++)
579  {
580  #define COVERAGE_IGNORE
581  if(verbose_during_solve)
582  {
583  std::cout << " Increment " << index << " of " << mpProblemDefinition->GetNumIncrementsForInitialDeformation() << "\n";
584  }
585  #undef COVERAGE_IGNORE
586 
587  if(mpProblemDefinition->GetTractionBoundaryConditionType()==PRESSURE_ON_DEFORMED)
588  {
589  mpProblemDefinition->SetPressureScaling(((double)index)/mpProblemDefinition->GetNumIncrementsForInitialDeformation());
590  }
591  mpMechanicsSolver->Solve();
592 
593  total_newton_iters += mpMechanicsSolver->GetNumNewtonIterations();
594  }
595 
596  mpMechanicsSolver->SetIncludeActiveTension(true);
597  MechanicsEventHandler::EndEvent(MechanicsEventHandler::ALL_MECH);
598  LOG(2, " Number of newton iterations = " << total_newton_iters);
599 
600 
601  unsigned mech_writer_counter = 0;
602 
603  if (mWriteOutput)
604  {
605  LOG(2, " Writing output");
606  mpMechanicsSolver->SetWriteOutput();
607  mpMechanicsSolver->WriteCurrentSpatialSolution("solution","nodes",mech_writer_counter);
608  p_cmgui_writer->WriteDeformationPositions(rGetDeformedPosition(), mech_writer_counter);
609 
610  if(!mNoElectricsOutput)
611  {
612  // the writer inside monodomain problem uses the printing timestep
613  // inside HeartConfig to estimate total number of timesteps, so make
614  // sure this is set to what we will use.
615  HeartConfig::Instance()->SetPrintingTimeStep(mpProblemDefinition->GetMechanicsSolveTimestep());
616 
617  // When we consider archiving these simulations we will need to get a bool back from the following
618  // command to decide whether or not the file is being extended, and hence whether to write the
619  // initial conditions into the .h5 file.
620  mpElectricsProblem->InitialiseWriter();
621  mpElectricsProblem->WriteOneStep(stepper.GetTime(), initial_voltage);
622  }
623 
624  if(mIsWatchedLocation)
625  {
626  WriteWatchedLocationData(stepper.GetTime(), initial_voltage);
627  }
628 
629  if(mNumTimestepsToOutputDeformationGradientsAndStress!=UNSIGNED_UNSET)
630  {
631  mpMechanicsSolver->WriteCurrentStrains(DEFORMATION_GRADIENT_F,"deformation_gradient",mech_writer_counter);
632  mpMechanicsSolver->WriteCurrentAverageElementStresses("second_PK",mech_writer_counter);
633  }
634  }
635 
636 
637  PrepareForSolve();
638 
642 // std::vector<double> current_solution_previous_time_step = mpMechanicsSolver->rGetCurrentSolution();
643 // std::vector<double> current_solution_second_last_time_step = mpMechanicsSolver->rGetCurrentSolution();
644 // bool first_step = true;
645 
646 
647  // reset this to false, may be reset again below
648  mpMechanicsSolver->SetComputeAverageStressPerElementDuringSolve(false);
649 
650  while (!stepper.IsTimeAtEnd())
651  {
652  LOG(2, "\nCurrent time = " << stepper.GetTime());
653  #define COVERAGE_IGNORE
654  if(verbose_during_solve)
655  {
656  // also output time to screen as newton solve information will be output
657  std::cout << "\n\n ** Current time = " << stepper.GetTime() << "\n";
658  }
659  #undef COVERAGE_IGNORE
660 
667  if(mpProblemDefinition->GetDeformationAffectsCellModels() || mpProblemDefinition->GetDeformationAffectsConductivity())
668  {
669  // Determine the stretch and deformation gradient on each element.
670  //
671  // If mpProblemDefinition->GetDeformationAffectsCellModels()==true:
672  // Stretch will be passed to the cell models.
673  //
674  // If mpProblemDefinition->GetDeformationAffectsConductivity()==true:
675  // The deformation gradient needs to be set up but does not need to be passed to the tissue
676  // so that F is used to compute the conductivity. Instead this is
677  // done through the line "mpElectricsProblem->GetMonodomainTissue()->SetConductivityModifier(this);" line above, which means
678  // rCalculateModifiedConductivityTensor() will be called on this class by the tissue, which then uses the F
679 
680  mpCardiacMechSolver->ComputeDeformationGradientAndStretchInEachElement(mDeformationGradientsForEachMechanicsElement, mStretchesForEachMechanicsElement);
681  }
682 
683  if( mpProblemDefinition->GetDeformationAffectsCellModels() )
684  {
685  // Set the stretches on each of the cell models
686  for(unsigned global_index = mpElectricsMesh->GetDistributedVectorFactory()->GetLow();
687  global_index < mpElectricsMesh->GetDistributedVectorFactory()->GetHigh();
688  global_index++)
689  {
690  unsigned containing_elem = mpMeshPair->rGetCoarseElementsForFineNodes()[global_index];
691  double stretch = mStretchesForEachMechanicsElement[containing_elem];
692  mpElectricsProblem->GetTissue()->GetCardiacCell(global_index)->SetStretch(stretch);
693  }
694  }
695 
696  p_electrics_solver->SetTimeStep(HeartConfig::Instance()->GetPdeTimeStep());
697 
703  LOG(2, " Solving electrics");
704  MechanicsEventHandler::BeginEvent(MechanicsEventHandler::NON_MECH);
705  for(unsigned i=0; i<mNumElecTimestepsPerMechTimestep; i++)
706  {
707  double current_time = stepper.GetTime() + i*HeartConfig::Instance()->GetPdeTimeStep();
708  double next_time = stepper.GetTime() + (i+1)*HeartConfig::Instance()->GetPdeTimeStep();
709 
710  // solve the electrics
711  p_electrics_solver->SetTimes(current_time, next_time);
712  p_electrics_solver->SetInitialCondition( initial_voltage );
713 
714  electrics_solution = p_electrics_solver->Solve();
715 
716  PetscReal min_voltage, max_voltage;
717  VecMax(electrics_solution,PETSC_NULL,&max_voltage); //the second param is where the index would be returned
718  VecMin(electrics_solution,PETSC_NULL,&min_voltage);
719  if(i==0)
720  {
721  LOG(2, " minimum and maximum voltage is " << min_voltage <<", "<<max_voltage);
722  }
723  else if(i==1)
724  {
725  LOG(2, " ..");
726  }
727 
728  PetscTools::Destroy(initial_voltage);
729  initial_voltage = electrics_solution;
730  }
731 
732  if(mpProblemDefinition->GetDeformationAffectsConductivity())
733  {
734  p_electrics_solver->SetMatrixIsNotAssembled();
735  }
736 
737 
743 
744  // compute Ca_I at each quad point (by interpolation, using the info on which
745  // electrics element the quad point is in. Then set Ca_I on the mechanics solver
746  LOG(2, " Interpolating Ca_I and voltage");
747 
748  //Collect the distributed calcium data into one Vec to be later replicated
749  for(unsigned node_index = 0; node_index<mpElectricsMesh->GetNumNodes(); node_index++)
750  {
751  if (mpElectricsMesh->GetDistributedVectorFactory()->IsGlobalIndexLocal(node_index))
752  {
753  double calcium_value = mpElectricsProblem->GetTissue()->GetCardiacCell(node_index)->GetIntracellularCalciumConcentration();
754  VecSetValue(calcium_data, node_index ,calcium_value, INSERT_VALUES);
755  }
756  }
757  PetscTools::Barrier();//not sure this is needed
758 
759  //Replicate electrics solution and calcium (replication is inside this constructor of ReplicatableVector)
760  ReplicatableVector electrics_solution_repl(electrics_solution);//size=(number of electrics nodes)*ELEC_PROB_DIM
761  ReplicatableVector calcium_repl(calcium_data);//size = number of electrics nodes
762 
763  //interpolate values onto mechanics mesh
764  for(unsigned i=0; i<mpMeshPair->rGetElementsAndWeights().size(); i++)
765  {
766  double interpolated_CaI = 0;
767  double interpolated_voltage = 0;
768 
769  Element<DIM,DIM>& element = *(mpElectricsMesh->GetElement(mpMeshPair->rGetElementsAndWeights()[i].ElementNum));
770 
771  for(unsigned node_index = 0; node_index<element.GetNumNodes(); node_index++)
772  {
773  unsigned global_index = element.GetNodeGlobalIndex(node_index);
774  double CaI_at_node = calcium_repl[global_index];
775  interpolated_CaI += CaI_at_node*mpMeshPair->rGetElementsAndWeights()[i].Weights(node_index);
776  //the following line assumes interleaved solution for ELEC_PROB_DIM>1 (e.g, [Vm_0, phi_e_0, Vm1, phi_e_1...])
777  interpolated_voltage += electrics_solution_repl[global_index*ELEC_PROB_DIM]*mpMeshPair->rGetElementsAndWeights()[i].Weights(node_index);
778  }
779 
780  assert(i<mInterpolatedCalciumConcs.size());
781  assert(i<mInterpolatedVoltages.size());
782  mInterpolatedCalciumConcs[i] = interpolated_CaI;
783  mInterpolatedVoltages[i] = interpolated_voltage;
784  }
785 
786  LOG(2, " Setting Ca_I. max value = " << Max(mInterpolatedCalciumConcs));
787 
788  // NOTE IF NHS: HERE WE SHOULD PERHAPS CHECK WHETHER THE CELL MODELS HAVE Ca_Trop
789  // AND UPDATE FROM NHS TO CELL_MODEL, BUT NOT SURE HOW TO DO THIS.. (esp for implicit)
790 
791  // set [Ca], V, t
792  mpCardiacMechSolver->SetCalciumAndVoltage(mInterpolatedCalciumConcs, mInterpolatedVoltages);
793  MechanicsEventHandler::EndEvent(MechanicsEventHandler::NON_MECH);
794 
795 
801  LOG(2, " Solving mechanics ");
802  mpMechanicsSolver->SetWriteOutput(false);
803 
804  // make sure the mechanics solver knows the current time (in case
805  // the traction say is time-dependent).
806  mpMechanicsSolver->SetCurrentTime(stepper.GetTime());
807 
808  // see if we will need to output stresses at the end of this timestep
809  if( mNumTimestepsToOutputDeformationGradientsAndStress!=UNSIGNED_UNSET
810  && (counter+1)%mNumTimestepsToOutputDeformationGradientsAndStress == 0 )
811  {
812  mpMechanicsSolver->SetComputeAverageStressPerElementDuringSolve(true);
813  }
814 
817 // for(unsigned i=0; i<mpMechanicsSolver->rGetCurrentSolution().size(); i++)
818 // {
819 // double current = mpMechanicsSolver->rGetCurrentSolution()[i];
820 // double previous = current_solution_previous_time_step[i];
821 // double second_last = current_solution_second_last_time_step[i];
822 // //double guess = 2*current - previous;
823 // double guess = 3*current - 3*previous + second_last;
824 //
825 // if(!first_step)
826 // {
827 // current_solution_second_last_time_step[i] = current_solution_previous_time_step[i];
828 // }
829 // current_solution_previous_time_step[i] = mpMechanicsSolver->rGetCurrentSolution()[i];
830 // mpMechanicsSolver->rGetCurrentSolution()[i] = guess;
831 // }
832 // first_step = false;
833 
834  MechanicsEventHandler::BeginEvent(MechanicsEventHandler::ALL_MECH);
835  mpCardiacMechSolver->Solve(stepper.GetTime(), stepper.GetNextTime(), mpProblemDefinition->GetContractionModelOdeTimestep());
836  MechanicsEventHandler::EndEvent(MechanicsEventHandler::ALL_MECH);
837 
838  LOG(2, " Number of newton iterations = " << mpMechanicsSolver->GetNumNewtonIterations());
839 
840  // update the current time
841  stepper.AdvanceOneTimeStep();
842  counter++;
843 
844 
850  MechanicsEventHandler::BeginEvent(MechanicsEventHandler::OUTPUT);
851  if(mWriteOutput && (counter%WRITE_EVERY_NTH_TIME==0))
852  {
853  LOG(2, " Writing output");
854  // write deformed position
855  mech_writer_counter++;
856  mpMechanicsSolver->SetWriteOutput();
857  mpMechanicsSolver->WriteCurrentSpatialSolution("solution","nodes",mech_writer_counter);
858 
859  p_cmgui_writer->WriteDeformationPositions(rGetDeformedPosition(), counter);
860 
861  if(!mNoElectricsOutput)
862  {
863  mpElectricsProblem->mpWriter->AdvanceAlongUnlimitedDimension();
864  mpElectricsProblem->WriteOneStep(stepper.GetTime(), electrics_solution);
865  }
866 
867  if(mIsWatchedLocation)
868  {
869  WriteWatchedLocationData(stepper.GetTime(), electrics_solution);
870  }
871  OnEndOfTimeStep(counter);
872 
873  if(mNumTimestepsToOutputDeformationGradientsAndStress!=UNSIGNED_UNSET && counter%mNumTimestepsToOutputDeformationGradientsAndStress==0)
874  {
875  mpMechanicsSolver->WriteCurrentStrains(DEFORMATION_GRADIENT_F,"deformation_gradient",mech_writer_counter);
876  mpMechanicsSolver->WriteCurrentAverageElementStresses("second_PK",mech_writer_counter);
877  }
878  mpMechanicsSolver->SetComputeAverageStressPerElementDuringSolve(false);
879  }
880  MechanicsEventHandler::EndEvent(MechanicsEventHandler::OUTPUT);
881 
882  // write the total elapsed time..
884  }
885 
886 
887 
888  if ((mWriteOutput) && (!mNoElectricsOutput))
889  {
891  mpElectricsProblem->mpWriter->Close();
892  delete mpElectricsProblem->mpWriter;
893 
894  std::string input_dir = mOutputDirectory+"/electrics";
895 
896  // Convert simulation data to meshalyzer format - commented
897  std::string config_directory = HeartConfig::Instance()->GetOutputDirectory();
899 
900 // Hdf5ToMeshalyzerConverter<DIM,DIM> meshalyzer_converter(FileFinder(input_dir, RelativeTo::ChasteTestOutput),
901 // "voltage", mpElectricsMesh,
902 // HeartConfig::Instance()->GetOutputUsingOriginalNodeOrdering(),
903 // HeartConfig::Instance()->GetVisualizerOutputPrecision() );
904 
905  // convert output to CMGUI format
907  "voltage", mpElectricsMesh, mHasBath,
908  HeartConfig::Instance()->GetVisualizerOutputPrecision());
909 
910  // Write mesh in a suitable form for meshalyzer
911  //std::string output_directory = mOutputDirectory + "/electrics/output";
912  // Write the mesh
913  //MeshalyzerMeshWriter<DIM,DIM> mesh_writer(output_directory, "mesh", false);
914  //mesh_writer.WriteFilesUsingMesh(*mpElectricsMesh);
915  // Write the parameters out
916  //HeartConfig::Instance()->Write();
917 
918  // interpolate the electrical data onto the mechanics mesh nodes and write CMGUI...
919  // Note: this calculates the data on ALL nodes of the mechanics mesh (incl internal,
920  // non-vertex ones), which won't be used if linear CMGUI visualisation
921  // of the mechanics solution is used.
922  VoltageInterpolaterOntoMechanicsMesh<DIM> converter(*mpElectricsMesh,*mpMechanicsMesh,variable_names,input_dir,"voltage");
923 
924  // reset to the default value
925  HeartConfig::Instance()->SetOutputDirectory(config_directory);
926  }
927 
928  if(p_cmgui_writer)
929  {
930  if(mNoElectricsOutput)
931  {
932  p_cmgui_writer->WriteCmguiScript("","undeformed");
933  }
934  else
935  {
936  p_cmgui_writer->WriteCmguiScript("../../electrics/cmgui_output/voltage_mechanics_mesh","undeformed");
937  }
938  delete p_cmgui_writer;
939  }
940  PetscTools::Destroy(electrics_solution);
941  PetscTools::Destroy(calcium_data);
942  delete p_electrics_solver;
943 
944  MechanicsEventHandler::EndEvent(MechanicsEventHandler::ALL);
945 }
946 
947 
948 
949 template<unsigned DIM, unsigned ELEC_PROB_DIM>
951 {
952  double max = -1e200;
953  for(unsigned i=0; i<vec.size(); i++)
954  {
955  if(vec[i]>max) max=vec[i];
956  }
957  return max;
958 }
959 
960 template<unsigned DIM, unsigned ELEC_PROB_DIM>
962 {
963  mNoElectricsOutput = true;
964 }
965 
966 template<unsigned DIM, unsigned ELEC_PROB_DIM>
968 {
969  mIsWatchedLocation = true;
970  mWatchedLocation = watchedLocation;
971 }
972 
973 template<unsigned DIM, unsigned ELEC_PROB_DIM>
975 {
976  mNumTimestepsToOutputDeformationGradientsAndStress = (unsigned) floor((timeStep/mpProblemDefinition->GetMechanicsSolveTimestep())+0.5);
977  if(fabs(mNumTimestepsToOutputDeformationGradientsAndStress*mpProblemDefinition->GetMechanicsSolveTimestep() - timeStep) > 1e-6)
978  {
979  EXCEPTION("Timestep provided for SetOutputDeformationGradientsAndStress() is not a multiple of mechanics solve timestep");
980  }
981 }
982 
983 
984 template<unsigned DIM, unsigned ELEC_PROB_DIM>
986 {
987  return mpMechanicsSolver->rGetDeformedPosition();
988 }
989 
990 
991 
992 
994 // Explicit instantiation
996 
997 //note: 1d incompressible material doesn't make sense
void Set(unsigned level, const std::string &rDirectory, const std::string &rFileName="log.txt")
Definition: LogFile.cpp:73
void SetOutputDirectory(const std::string &rOutputDirectory)
void WriteWatchedLocationData(double time, Vec voltage)
static void Barrier(const std::string callerId="")
Definition: PetscTools.cpp:134
boost::numeric::ublas::c_matrix< T, 1, 1 > Inverse(const boost::numeric::ublas::c_matrix< T, 1, 1 > &rM)
unsigned GetNodeGlobalIndex(unsigned localIndex) const
#define EXCEPTION(message)
Definition: Exception.hpp:143
double Max(std::vector< double > &vec)
static AbstractCardiacProblem< DIM, DIM, 2u > * Create(ElectricsProblemType problemType, AbstractCardiacCellFactory< DIM > *pCellFactory)
double GetPdeTimeStep() const
#define NEVER_REACHED
Definition: Exception.hpp:206
bool OptionExists(const std::string &rOption)
void SetTimes(double tStart, double tEnd)
std::vector< c_vector< double, DIM > > & rGetDeformedPosition()
out_stream OpenOutputFile(const std::string &rFileName, std::ios_base::openmode mode=std::ios::out|std::ios::trunc) const
CardiacElectroMechanicsProblem(CompressibilityType compressibilityType, ElectricsProblemType electricsProblemType, TetrahedralMesh< DIM, DIM > *pElectricsMesh, QuadraticMesh< DIM > *pMechanicsMesh, AbstractCardiacCellFactory< DIM > *pCellFactory, ElectroMechanicsProblemDefinition< DIM > *pProblemDefinition, std::string outputDirectory)
static void Close()
Definition: LogFile.cpp:101
unsigned GetNumNodes() const
const unsigned UNSIGNED_UNSET
Definition: Exception.hpp:52
void SetOutputDeformationGradientsAndStress(double timestep)
static AbstractCardiacProblem< DIM, DIM, 1u > * Create(ElectricsProblemType problemType, AbstractCardiacCellFactory< DIM > *pCellFactory)
static void Destroy(Vec &rVec)
Definition: PetscTools.hpp:351
c_matrix< double, DIM, DIM > & rCalculateModifiedConductivityTensor(unsigned elementIndex, const c_matrix< double, DIM, DIM > &rOriginalConductivity, unsigned domainIndex)
void WriteHeader(std::string simulationType="")
Definition: LogFile.cpp:111
AbstractCardiacProblem< DIM, DIM, ELEC_PROB_DIM > * mpElectricsProblem
void WriteElapsedTime(std::string pre="")
Definition: LogFile.cpp:116
static CommandLineArguments * Instance()
void SetUpBoxesOnFineMesh(double boxWidth=-1)
std::string GetOutputDirectory() const
static LogFile * Instance()
Definition: LogFile.cpp:52
virtual void WriteFilesUsingMesh(AbstractTetrahedralMesh< ELEMENT_DIM, SPACE_DIM > &rMesh, bool keepOriginalElementIndexing=true)
static void Reset()
void SetPrintingTimeStep(double printingTimeStep)
void SetOutputFilenamePrefix(const std::string &rOutputFilenamePrefix)
static HeartConfig * Instance()
static AbstractCardiacProblem< DIM, DIM, PROBLEM_DIM > * Create(ElectricsProblemType problemType, AbstractCardiacCellFactory< DIM > *pCellFactory)
void SetWatchedPosition(c_vector< double, DIM > watchedLocation)