347 assert(mpElectricsMesh!=NULL);
348 assert(mpMechanicsMesh!=NULL);
349 assert(mpProblemDefinition!=NULL);
350 assert(mpCardiacMechSolver==NULL);
355 if (fabs(mNumElecTimestepsPerMechTimestep*
HeartConfig::Instance()->GetPdeTimeStep() - mpProblemDefinition->GetMechanicsSolveTimestep()) > 1e-6)
357 EXCEPTION(
"Electrics PDE timestep does not divide mechanics solve timestep");
362 std::string log_dir = mOutputDirectory;
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]");
370 LOG(2,
"Electrics mesh has " << mpElectricsMesh->GetNumNodes() <<
" nodes");
371 LOG(2,
"Mechanics mesh has " << mpMechanicsMesh->GetNumNodes() <<
" nodes");
373 LOG(2,
"Initialising..");
376 if (mIsWatchedLocation)
378 DetermineWatchedNodes();
382 mpElectricsProblem->SetMesh(mpElectricsMesh);
383 mpElectricsProblem->Initialise();
385 if (mCompressibilityType==INCOMPRESSIBLE)
387 switch(mpProblemDefinition->GetSolverType())
391 *mpMechanicsMesh,*mpProblemDefinition,mDeformationOutputDirectory);
395 *mpMechanicsMesh,*mpProblemDefinition,mDeformationOutputDirectory);
405 switch(mpProblemDefinition->GetSolverType())
409 *mpMechanicsMesh,*mpProblemDefinition,mDeformationOutputDirectory);
413 *mpMechanicsMesh,*mpProblemDefinition,mDeformationOutputDirectory);
422 assert(mpMechanicsSolver);
427 mpMeshPair->SetUpBoxesOnFineMesh();
428 mpMeshPair->ComputeFineElementsAndWeightsForCoarseQuadPoints(*(mpCardiacMechSolver->GetQuadratureRule()),
false);
429 mpMeshPair->DeleteFineBoxCollection();
431 mpCardiacMechSolver->SetFineCoarseMeshPair(mpMeshPair);
432 mpCardiacMechSolver->Initialise();
434 unsigned num_quad_points = mpCardiacMechSolver->GetTotalNumQuadPoints();
435 mInterpolatedCalciumConcs.assign(num_quad_points, 0.0);
436 mInterpolatedVoltages.assign(num_quad_points, 0.0);
438 if (mpProblemDefinition->ReadFibreSheetDirectionsFromFile())
440 mpCardiacMechSolver->SetVariableFibreSheetDirections(mpProblemDefinition->GetFibreSheetDirectionsFile(),
441 mpProblemDefinition->GetFibreSheetDirectionsDefinedPerQuadraturePoint());
445 if (mpProblemDefinition->GetDeformationAffectsConductivity() || mpProblemDefinition->GetDeformationAffectsCellModels())
447 mpMeshPair->SetUpBoxesOnCoarseMesh();
451 if (mpProblemDefinition->GetDeformationAffectsCellModels() || mpProblemDefinition->GetDeformationAffectsConductivity())
454 mStretchesForEachMechanicsElement.resize(mpMechanicsMesh->GetNumElements(), 1.0);
457 mDeformationGradientsForEachMechanicsElement.resize(mpMechanicsMesh->GetNumElements(),identity_matrix<double>(DIM));
461 if (mpProblemDefinition->GetDeformationAffectsCellModels())
465 mpMeshPair->ComputeCoarseElementsForFineNodes(
false);
469 if (mpProblemDefinition->GetDeformationAffectsConductivity())
473 mpMeshPair->ComputeCoarseElementsForFineElementCentroids(
false);
477 mpElectricsProblem->GetTissue()->SetConductivityModifier(
this);
491 if (mpCardiacMechSolver==NULL)
496 bool verbose_during_solve = ( mpProblemDefinition->GetVerboseDuringSolve()
501 mpProblemDefinition->Validate();
504 p_bcc->DefineZeroNeumannOnMeshBoundary(mpElectricsMesh, 0);
505 mpElectricsProblem->SetBoundaryConditionsContainer(p_bcc);
510 = mpElectricsProblem->CreateSolver();
513 Vec electrics_solution=NULL;
514 Vec calcium_data= mpElectricsMesh->GetDistributedVectorFactory()->CreateVec();
515 Vec initial_voltage = mpElectricsProblem->CreateInitialCondition();
518 unsigned counter = 0;
524 std::vector<std::string> variable_names;
528 mpMechanicsSolver->SetWriteOutput();
529 mpMechanicsSolver->WriteCurrentSpatialSolution(
"undeformed",
"nodes");
533 *(this->mpMechanicsMesh),
534 WRITE_QUADRATIC_MESH);
535 variable_names.push_back(
"V");
536 if (ELEC_PROB_DIM==2)
538 variable_names.push_back(
"Phi_e");
541 std::vector<std::string> regions;
542 regions.push_back(
"tissue");
543 regions.push_back(
"bath");
558 LOG(2,
"\nSolving for initial deformation");
560 if (verbose_during_solve)
562 std::cout <<
"\n\n ** Solving for initial deformation\n";
566 mpMechanicsSolver->SetWriteOutput(
false);
568 mpMechanicsSolver->SetCurrentTime(0.0);
571 mpMechanicsSolver->SetIncludeActiveTension(
false);
572 if (mNumTimestepsToOutputDeformationGradientsAndStress!=
UNSIGNED_UNSET)
574 mpMechanicsSolver->SetComputeAverageStressPerElementDuringSolve(
true);
577 unsigned total_newton_iters = 0;
578 for (
unsigned index=1; index<=mpProblemDefinition->GetNumIncrementsForInitialDeformation(); index++)
581 if (verbose_during_solve)
583 std::cout <<
" Increment " << index <<
" of " << mpProblemDefinition->GetNumIncrementsForInitialDeformation() <<
"\n";
587 if (mpProblemDefinition->GetTractionBoundaryConditionType()==PRESSURE_ON_DEFORMED)
589 mpProblemDefinition->SetPressureScaling(((
double)index)/mpProblemDefinition->GetNumIncrementsForInitialDeformation());
591 mpMechanicsSolver->Solve();
593 total_newton_iters += mpMechanicsSolver->GetNumNewtonIterations();
596 mpMechanicsSolver->SetIncludeActiveTension(
true);
598 LOG(2,
" Number of newton iterations = " << total_newton_iters);
601 unsigned mech_writer_counter = 0;
605 LOG(2,
" Writing output");
606 mpMechanicsSolver->SetWriteOutput();
607 mpMechanicsSolver->WriteCurrentSpatialSolution(
"solution",
"nodes",mech_writer_counter);
610 if (!mNoElectricsOutput)
620 mpElectricsProblem->InitialiseWriter();
621 mpElectricsProblem->WriteOneStep(stepper.
GetTime(), initial_voltage);
624 if (mIsWatchedLocation)
626 WriteWatchedLocationData(stepper.
GetTime(), initial_voltage);
629 if (mNumTimestepsToOutputDeformationGradientsAndStress!=
UNSIGNED_UNSET)
631 mpMechanicsSolver->WriteCurrentStrains(DEFORMATION_GRADIENT_F,
"deformation_gradient",mech_writer_counter);
632 mpMechanicsSolver->WriteCurrentAverageElementStresses(
"second_PK",mech_writer_counter);
648 mpMechanicsSolver->SetComputeAverageStressPerElementDuringSolve(
false);
652 LOG(2,
"\nCurrent time = " << stepper.
GetTime());
654 if (verbose_during_solve)
657 std::cout <<
"\n\n ** Current time = " << stepper.
GetTime() <<
"\n";
667 if (mpProblemDefinition->GetDeformationAffectsCellModels() || mpProblemDefinition->GetDeformationAffectsConductivity())
680 mpCardiacMechSolver->ComputeDeformationGradientAndStretchInEachElement(mDeformationGradientsForEachMechanicsElement, mStretchesForEachMechanicsElement);
683 if (mpProblemDefinition->GetDeformationAffectsCellModels())
686 for (
unsigned global_index = mpElectricsMesh->GetDistributedVectorFactory()->GetLow();
687 global_index < mpElectricsMesh->GetDistributedVectorFactory()->GetHigh();
690 unsigned containing_elem = mpMeshPair->rGetCoarseElementsForFineNodes()[global_index];
691 double stretch = mStretchesForEachMechanicsElement[containing_elem];
692 mpElectricsProblem->GetTissue()->GetCardiacCell(global_index)->SetStretch(stretch);
703 LOG(2,
" Solving electrics");
705 for (
unsigned i=0; i<mNumElecTimestepsPerMechTimestep; i++)
711 p_electrics_solver->
SetTimes(current_time, next_time);
714 electrics_solution = p_electrics_solver->
Solve();
716 PetscReal min_voltage, max_voltage;
721 LOG(2,
" minimum and maximum voltage is " << min_voltage <<
", "<<max_voltage);
729 initial_voltage = electrics_solution;
732 if (mpProblemDefinition->GetDeformationAffectsConductivity())
746 LOG(2,
" Interpolating Ca_I and voltage");
749 for (
unsigned node_index = 0; node_index<mpElectricsMesh->GetNumNodes(); node_index++)
751 if (mpElectricsMesh->GetDistributedVectorFactory()->IsGlobalIndexLocal(node_index))
753 double calcium_value = mpElectricsProblem->GetTissue()->GetCardiacCell(node_index)->GetIntracellularCalciumConcentration();
754 VecSetValue(calcium_data, node_index ,calcium_value, INSERT_VALUES);
764 for (
unsigned i=0; i<mpMeshPair->rGetElementsAndWeights().size(); i++)
766 double interpolated_CaI = 0;
767 double interpolated_voltage = 0;
769 Element<DIM,DIM>& element = *(mpElectricsMesh->GetElement(mpMeshPair->rGetElementsAndWeights()[i].ElementNum));
771 for (
unsigned node_index = 0; node_index<element.
GetNumNodes(); node_index++)
774 double CaI_at_node = calcium_repl[global_index];
775 interpolated_CaI += CaI_at_node*mpMeshPair->rGetElementsAndWeights()[i].Weights(node_index);
777 interpolated_voltage += electrics_solution_repl[global_index*ELEC_PROB_DIM]*mpMeshPair->rGetElementsAndWeights()[i].Weights(node_index);
780 assert(i<mInterpolatedCalciumConcs.size());
781 assert(i<mInterpolatedVoltages.size());
782 mInterpolatedCalciumConcs[i] = interpolated_CaI;
783 mInterpolatedVoltages[i] = interpolated_voltage;
786 LOG(2,
" Setting Ca_I. max value = " << Max(mInterpolatedCalciumConcs));
792 mpCardiacMechSolver->SetCalciumAndVoltage(mInterpolatedCalciumConcs, mInterpolatedVoltages);
801 LOG(2,
" Solving mechanics ");
802 mpMechanicsSolver->SetWriteOutput(
false);
806 mpMechanicsSolver->SetCurrentTime(stepper.
GetTime());
809 if (mNumTimestepsToOutputDeformationGradientsAndStress!=
UNSIGNED_UNSET
810 && (counter+1)%mNumTimestepsToOutputDeformationGradientsAndStress == 0)
812 mpMechanicsSolver->SetComputeAverageStressPerElementDuringSolve(
true);
835 mpCardiacMechSolver->Solve(stepper.
GetTime(), stepper.
GetNextTime(), mpProblemDefinition->GetContractionModelOdeTimestep());
838 LOG(2,
" Number of newton iterations = " << mpMechanicsSolver->GetNumNewtonIterations());
851 if (mWriteOutput && (counter%WRITE_EVERY_NTH_TIME==0))
853 LOG(2,
" Writing output");
855 mech_writer_counter++;
856 mpMechanicsSolver->SetWriteOutput();
857 mpMechanicsSolver->WriteCurrentSpatialSolution(
"solution",
"nodes",mech_writer_counter);
861 if (!mNoElectricsOutput)
863 mpElectricsProblem->mpWriter->AdvanceAlongUnlimitedDimension();
864 mpElectricsProblem->WriteOneStep(stepper.
GetTime(), electrics_solution);
867 if (mIsWatchedLocation)
869 WriteWatchedLocationData(stepper.
GetTime(), electrics_solution);
871 OnEndOfTimeStep(counter);
873 if (mNumTimestepsToOutputDeformationGradientsAndStress!=
UNSIGNED_UNSET && counter%mNumTimestepsToOutputDeformationGradientsAndStress==0)
875 mpMechanicsSolver->WriteCurrentStrains(DEFORMATION_GRADIENT_F,
"deformation_gradient",mech_writer_counter);
876 mpMechanicsSolver->WriteCurrentAverageElementStresses(
"second_PK",mech_writer_counter);
878 mpMechanicsSolver->SetComputeAverageStressPerElementDuringSolve(
false);
886 if ((mWriteOutput) && (!mNoElectricsOutput))
889 mpElectricsProblem->mpWriter->Close();
890 delete mpElectricsProblem->mpWriter;
892 std::string input_dir = mOutputDirectory+
"/electrics";
905 "voltage", mpElectricsMesh, mHasBath,
928 if (mNoElectricsOutput)
934 p_cmgui_writer->
WriteCmguiScript(
"../../electrics/cmgui_output/voltage_mechanics_mesh",
"undeformed");
936 delete p_cmgui_writer;
940 delete p_electrics_solver;