Chaste Commit::baa90ac2819b962188b7562f2326be23c47859a7
CardiacElectroMechanicsProblem.cpp
1/*
2
3Copyright (c) 2005-2024, University of Oxford.
4All rights reserved.
5
6University of Oxford means the Chancellor, Masters and Scholars of the
7University of Oxford, having an administrative office at Wellington
8Square, Oxford OX1 2JD, UK.
9
10This file is part of Chaste.
11
12Redistribution and use in source and binary forms, with or without
13modification, 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
23THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
24AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
25IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
26ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
27LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
28CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE
29GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
30HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
31LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
32OF 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
59template<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 // LCOV_EXCL_START
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 // LCOV_EXCL_STOP
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 // LCOV_EXCL_START
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 // LCOV_EXCL_STOP
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
143template<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
168template<unsigned DIM, unsigned ELEC_PROB_DIM>
169c_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
188template<unsigned DIM, unsigned PROBLEM_DIM>
190{
191public:
198 static AbstractCardiacProblem<DIM, DIM, PROBLEM_DIM>* Create(ElectricsProblemType problemType,
199 AbstractCardiacCellFactory<DIM>* pCellFactory);
200};
201
205template<unsigned DIM>
207{
208public:
215 static AbstractCardiacProblem<DIM, DIM, 1u>* Create(ElectricsProblemType problemType,
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
229template<unsigned DIM>
231{
232public:
239 static AbstractCardiacProblem<DIM, DIM, 2u>* Create(ElectricsProblemType problemType,
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
255template<unsigned DIM, unsigned ELEC_PROB_DIM>
257 CompressibilityType compressibilityType,
258 ElectricsProblemType electricsProblemType,
259 TetrahedralMesh<DIM,DIM>* pElectricsMesh,
260 QuadraticMesh<DIM>* pMechanicsMesh,
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
326template<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
342}
343
344template<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
487template<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 // LCOV_EXCL_START
560 if (verbose_during_solve)
561 {
562 std::cout << "\n\n ** Solving for initial deformation\n";
563 }
564 // LCOV_EXCL_STOP
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 // LCOV_EXCL_START
581 if (verbose_during_solve)
582 {
583 std::cout << " Increment " << index << " of " << mpProblemDefinition->GetNumIncrementsForInitialDeformation() << "\n";
584 }
585 // LCOV_EXCL_STOP
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 // LCOV_EXCL_START
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 // LCOV_EXCL_STOP
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,CHASTE_PETSC_NULLPTR,&max_voltage); //the second param is where the index would be returned
718 VecMin(electrics_solution,CHASTE_PETSC_NULLPTR,&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 if ((mWriteOutput) && (!mNoElectricsOutput))
887 {
889 mpElectricsProblem->mpWriter->Close();
890 delete mpElectricsProblem->mpWriter;
891
892 std::string input_dir = mOutputDirectory+"/electrics";
893
894 // Convert simulation data to meshalyzer format - commented
895 std::string config_directory = HeartConfig::Instance()->GetOutputDirectory();
897
898// Hdf5ToMeshalyzerConverter<DIM,DIM> meshalyzer_converter(FileFinder(input_dir, RelativeTo::ChasteTestOutput),
899// "voltage", mpElectricsMesh,
900// HeartConfig::Instance()->GetOutputUsingOriginalNodeOrdering(),
901// HeartConfig::Instance()->GetVisualizerOutputPrecision() );
902
903 // convert output to CMGUI format
905 "voltage", mpElectricsMesh, mHasBath,
906 HeartConfig::Instance()->GetVisualizerOutputPrecision());
907
908 // Write mesh in a suitable form for meshalyzer
909 //std::string output_directory = mOutputDirectory + "/electrics/output";
910 // Write the mesh
911 //MeshalyzerMeshWriter<DIM,DIM> mesh_writer(output_directory, "mesh", false);
912 //mesh_writer.WriteFilesUsingMesh(*mpElectricsMesh);
913 // Write the parameters out
914 //HeartConfig::Instance()->Write();
915
916 // interpolate the electrical data onto the mechanics mesh nodes and write CMGUI...
917 // Note: this calculates the data on ALL nodes of the mechanics mesh (incl internal,
918 // non-vertex ones), which won't be used if linear CMGUI visualisation
919 // of the mechanics solution is used.
920 VoltageInterpolaterOntoMechanicsMesh<DIM> converter(*mpElectricsMesh,*mpMechanicsMesh,variable_names,input_dir,"voltage");
921
922 // reset to the default value
923 HeartConfig::Instance()->SetOutputDirectory(config_directory);
924 }
925
926 if (p_cmgui_writer)
927 {
928 if (mNoElectricsOutput)
929 {
930 p_cmgui_writer->WriteCmguiScript("","undeformed");
931 }
932 else
933 {
934 p_cmgui_writer->WriteCmguiScript("../../electrics/cmgui_output/voltage_mechanics_mesh","undeformed");
935 }
936 delete p_cmgui_writer;
937 }
938 PetscTools::Destroy(electrics_solution);
939 PetscTools::Destroy(calcium_data);
940 delete p_electrics_solver;
941
942 MechanicsEventHandler::EndEvent(MechanicsEventHandler::ALL);
943}
944
945template<unsigned DIM, unsigned ELEC_PROB_DIM>
947{
948 double max = -1e200;
949 for (unsigned i=0; i<vec.size(); i++)
950 {
951 if (vec[i]>max) max=vec[i];
952 }
953 return max;
954}
955
956template<unsigned DIM, unsigned ELEC_PROB_DIM>
961
962template<unsigned DIM, unsigned ELEC_PROB_DIM>
964{
965 mIsWatchedLocation = true;
966 mWatchedLocation = watchedLocation;
967}
968
969template<unsigned DIM, unsigned ELEC_PROB_DIM>
971{
972 mNumTimestepsToOutputDeformationGradientsAndStress = (unsigned) floor((timeStep/mpProblemDefinition->GetMechanicsSolveTimestep())+0.5);
973 if (fabs(mNumTimestepsToOutputDeformationGradientsAndStress*mpProblemDefinition->GetMechanicsSolveTimestep() - timeStep) > 1e-6)
974 {
975 EXCEPTION("Timestep provided for SetOutputDeformationGradientsAndStress() is not a multiple of mechanics solve timestep");
976 }
977}
978
979template<unsigned DIM, unsigned ELEC_PROB_DIM>
981{
982 return mpMechanicsSolver->rGetDeformedPosition();
983}
984
985// Explicit instantiation
986
987//note: 1d incompressible material doesn't make sense
#define EXCEPTION(message)
#define NEVER_REACHED
const unsigned UNSIGNED_UNSET
Definition Exception.hpp:53
#define CHASTE_PETSC_NULLPTR
A macro to define PETSc null pointer based on the PETSc version.
boost::numeric::ublas::c_matrix< T, 1, 1 > Inverse(const boost::numeric::ublas::c_matrix< T, 1, 1 > &rM)
void SetTimes(double tStart, double tEnd)
unsigned GetNumNodes() const
unsigned GetNodeGlobalIndex(unsigned localIndex) const
virtual void WriteFilesUsingMesh(AbstractTetrahedralMesh< ELEMENT_DIM, SPACE_DIM > &rMesh, bool keepOriginalElementIndexing=true)
double Max(std::vector< double > &vec)
CardiacElectroMechanicsProblem(CompressibilityType compressibilityType, ElectricsProblemType electricsProblemType, TetrahedralMesh< DIM, DIM > *pElectricsMesh, QuadraticMesh< DIM > *pMechanicsMesh, AbstractCardiacCellFactory< DIM > *pCellFactory, ElectroMechanicsProblemDefinition< DIM > *pProblemDefinition, std::string outputDirectory)
c_matrix< double, DIM, DIM > & rCalculateModifiedConductivityTensor(unsigned elementIndex, const c_matrix< double, DIM, DIM > &rOriginalConductivity, unsigned domainIndex)
void SetWatchedPosition(c_vector< double, DIM > watchedLocation)
std::vector< c_vector< double, DIM > > & rGetDeformedPosition()
void WriteWatchedLocationData(double time, Vec voltage)
AbstractCardiacProblem< DIM, DIM, ELEC_PROB_DIM > * mpElectricsProblem
void WriteCmguiScript(std::string fieldBaseName="", std::string undeformedBaseName="")
void WriteDeformationPositions(std::vector< c_vector< double, DIM > > &rDeformedPositions, unsigned counter)
void WriteInitialMesh(std::string fileName="")
void SetAdditionalFieldNames(std::vector< std::string > &rFieldNames)
void SetRegionNames(std::vector< std::string > &rRegionNames)
bool OptionExists(const std::string &rOption)
static CommandLineArguments * Instance()
static AbstractCardiacProblem< DIM, DIM, 1u > * Create(ElectricsProblemType problemType, AbstractCardiacCellFactory< DIM > *pCellFactory)
static AbstractCardiacProblem< DIM, DIM, 2u > * Create(ElectricsProblemType problemType, AbstractCardiacCellFactory< DIM > *pCellFactory)
static AbstractCardiacProblem< DIM, DIM, PROBLEM_DIM > * Create(ElectricsProblemType problemType, AbstractCardiacCellFactory< DIM > *pCellFactory)
double GetPdeTimeStep() const
void SetPrintingTimeStep(double printingTimeStep)
void SetOutputDirectory(const std::string &rOutputDirectory)
void SetOutputFilenamePrefix(const std::string &rOutputFilenamePrefix)
static void Reset()
std::string GetOutputDirectory() const
static HeartConfig * Instance()
void Set(unsigned level, const std::string &rDirectory, const std::string &rFileName="log.txt")
Definition LogFile.cpp:73
static LogFile * Instance()
Definition LogFile.cpp:52
void WriteElapsedTime(std::string pre="")
Definition LogFile.cpp:116
void WriteHeader(std::string simulationType="")
Definition LogFile.cpp:111
static void Close()
Definition LogFile.cpp:101
out_stream OpenOutputFile(const std::string &rFileName, std::ios_base::openmode mode=std::ios::out|std::ios::trunc) const
static void Destroy(Vec &rVec)
static void Barrier(const std::string callerId="")
bool IsTimeAtEnd() const
double GetTime() const
void AdvanceOneTimeStep()
double GetNextTime() const