PacemakerBidomainSolver.cpp

Go to the documentation of this file.

/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
   This file is part of CEPS.
 
   CEPS is free software: you can redistribute it and/or modify
   it under the terms of the GNU General Public License as published by
   the Free Software Foundation, either version 3 of the License, or
   (at your option) any later version.
 
   CEPS is distributed in the hope that it will be useful,
   but WITHOUT ANY WARRANTY; without even the implied warranty of
   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
   GNU General Public License for more details.
 
   You should have received a copy of the GNU General Public License
   along with CEPS (see file LICENSE at root of project).
   If not, see <https://www.gnu.org/licenses/>.
 
 
   Copyright 2019-2024 Inria, Universite de Bordeaux
 
   Authors, in alphabetical order:
 
     Pierre-Elliott BECUE, Florian CARO, Yves COUDIERE(*), Andjela DAVIDOVIC, 
     Charlie DOUANLA-LONTSI, Marc FUENTES, Mehdi JUHOOR, Michael LEGUEBE(*), 
     Pauline MIGERDITICHAN, Valentin PANNETIER(*), Nejib ZEMZEMI.
     * : currently active authors
 
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
#include "cardiac/pacemaker/solver/PacemakerBidomainSolver.hpp"
 
#include "cardiac/assembler/FE/FECardiacBCAssembler.hpp"
#include "common/toolbox/VariableTimeStepper.hpp"
#include "cardiac/pacemaker/assembler/FE/FEPacemakerBidomainAssembler.hpp"
#include "cardiac/pacemaker/problem/PacemakerBidomainProblem.hpp"
#include "pde/assemblers/FE/FEIntegrator.hpp"
 
PacemakerBidomainSolver::PacemakerBidomainSolver(PacemakerBidomainProblem* problem) :
    ExtendedBidomainSolver (problem),
    AbstractPacemakerSolver(problem,problem->getParameters())
{
  // renew the time stepper with a variable one
  ceps::destroyObject(m_timeStepper);
  m_timeStepper = ceps::getNew<VariableTimeStepper>(
    problem->getStartTime(), problem->getEndTime(), problem->getTimeStep()
  );
 
  if (ceps::isValidPtr(m_timeWriter))
  {
    m_timeWriter->addCustomData("Current I (uA)", &m_current);
    m_timeWriter->addCustomData("Measured voltage (mV)", &m_Vcroco);
    m_timeWriter->addCustomData("Mean of ue on anode (mV)", &m_meanOfPotentialOnAnode);
    m_timeWriter->addCustomData("Mean of ue on cathode (mV)", &m_meanOfPotentialOnCathode);
    m_timeWriter->addCustomData("Phase", &m_phase);
  }
}
 
void
PacemakerBidomainSolver::initializeAssemblers()
{
  auto pbpb = getPacemakerBidomainProblem();
  m_opAsb = ceps::getNew<FEPacemakerBidomainAssembler>(pbpb, m_fe);
  m_bcAsb = ceps::getNew<FECardiacBCAssembler>(pbpb,CepsSet<CepsAttribute>({CepsUniversal}),m_fe);
}
 
PacemakerBidomainProblem*
PacemakerBidomainSolver::getPacemakerBidomainProblem() const
{
  return ceps::runtimeCast<PacemakerBidomainProblem*>(getTimedProblem());
}
 
void
PacemakerBidomainSolver::assembleAndSolve()
{
 
  CepsReal time     = m_timeStepper->getTime();
  CepsReal previous = m_timeStepper->getTimeStep();
  CepsReal dt       = getRequiredTimeStep(time,getPacemakerBidomainProblem()->getTimeStep());
 
  if (not ceps::equals(previous,dt))
    m_timeStepper->setTimeStep(dt);
 
  ExtendedBidomainSolver::assembleAndSolve();
 
  // Compute addition values
  auto pb                    = getPacemakerBidomainProblem();
  auto anode                 = pb->getAnode();
  auto cathode               = pb->getCathode();
  auto ue                    = pb->getUeUnknown();
  m_current                  = getPacemakerCurrent(m_uNp1, m_fe, time + dt);
  m_Vcroco                   = getAlligatorVoltage(m_uNp1, m_fe, time + dt);
  m_meanOfPotentialOnAnode   = getMeanOnElectrodeOf(m_uNp1, ue, m_fe, anode);
  m_meanOfPotentialOnCathode = getMeanOnElectrodeOf(m_uNp1, ue, m_fe, cathode);
  m_phase                    = static_cast<CepsReal>(pb->getPhase(time));
  #ifdef CEPS_DEBUG_ENABLED
    ceps::checkNanOrInf(m_current                 ,"pacemaker current"     );
    ceps::checkNanOrInf(m_Vcroco                  ,"voltage between clips" );
    ceps::checkNanOrInf(m_meanOfPotentialOnAnode  ,"mean anode potential"  );
    ceps::checkNanOrInf(m_meanOfPotentialOnCathode,"mean cathode potential");
  #endif
}