PacemakerBidomainProblem.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/problem/PacemakerBidomainProblem.hpp"
 
#include "cardiac/pacemaker/solver/PacemakerBidomainSolver.hpp"
#include "pde/problem/static/FluxAnodeCathodeProblem.hpp"
#include "pde/solver/static/FluxAnodeCathodeSolver.hpp"
 
PacemakerBidomainProblem::PowerSupply::PowerSupply(PacemakerBidomainProblem* problem) :
    AbstractPacemakerProblem::PowerSupply(problem),
    m_problem(problem)
{}
 
CepsReal
PacemakerBidomainProblem::PowerSupply::eval(CepsStandardArgs args)
{
  PacemakerBidomainProblem* pb = ceps::runtimeCast<PacemakerBidomainProblem*>(m_problem);
 
  if (args.unknownId == pb->getVCetsUnknown()->getIdentifier())
  {
    auto     Am  = pb->getAm()->getFunctor()->eval(args);
    auto     Cm  = pb->getCm()->getFunctor()->eval(args);
    CepsReal res = AbstractPacemakerProblem::PowerSupply::eval(args);
 
    return res / (Am * Cm);
  }
  return 0.0;
}
 
PacemakerBidomainProblem::PacemakerBidomainProblem(Geometry* g, InputParameters* params) :
    ExtendedBidomainProblem (g, params),
    AbstractPacemakerProblem(g, params)
{
  setProblemName("Extended Bidomain connected to a Pacemaker");
 
  if (ceps::isValidPtr(params))
    PacemakerBidomainProblem::setupWithParameters(params);
}
 
void
PacemakerBidomainProblem::getInitialCondition(DHVecPtr v0) const
{
  v0->zero();
  ExtendedBidomainProblem::getInitialCondition(v0);
 
  // auto QCetsDof = m_discr->getDegreesOfFreedomForUnknown(m_unknowns[4])[0];
  // if (QCetsDof->getOwner() == ceps::getRank())
  // {
  //   v0->getLocalData();
  //   (*v0)[QCetsDof->getGlobalIndex()] = m_pmkCets * m_stimAmplitude;
  //   v0->releaseLocalData();
  // }
}
 
void
PacemakerBidomainProblem::defineUnknowns()
{
  addUnknown(UnknownsName::vm, getTissueAttributes     ());
  addUnknown(UnknownsName::ue, getWholeDomainAttributes());  // ue is defined everywhere
  addUnknownInteraction(UnknownsName::vm, UnknownsName::ue, getTissueAttributes());
 
  // the 0D unknowns "potential"
  addZeroDUnknown(UnknownsName::UAnode  );
  addZeroDUnknown(UnknownsName::UCathode);
  addUnknownInteraction(UnknownsName::UAnode, UnknownsName::UCathode);
 
  // and the 0D unknowns
  addZeroDUnknown(UnknownsName::VCets);
  addZeroDUnknown(UnknownsName::VCts );
  addUnknownInteraction(UnknownsName::VCets   , UnknownsName::VCts    );
  addUnknownInteraction(UnknownsName::VCets   , UnknownsName::UAnode  );
  addUnknownInteraction(UnknownsName::VCets   , UnknownsName::UCathode);
  addUnknownInteraction(UnknownsName::VCts    , UnknownsName::UAnode  );
  addUnknownInteraction(UnknownsName::VCts    , UnknownsName::UCathode);
  // sigma_e U_e.gradN is integrated on the cathode to get the current
  addUnknownInteraction(UnknownsName::UCathode, UnknownsName::ue, getCathode()->getAttributes());
  addUnknownInteraction(UnknownsName::UAnode  , UnknownsName::ue, getCathode()->getAttributes());
}
 
void
PacemakerBidomainProblem::defineBoundaryConditions()
{}
 
void
PacemakerBidomainProblem::defineSourceTerms()
{
  ExtendedBidomainProblem::defineSourceTerms();
 
  m_functions->addEntry("capaChargeFunc", ceps::getNew<PowerSupply>(this));
  m_sourceTerms->add(
    "__power_supply capaChargeFunc UNKNOWN "  //
    + ceps::toString(getVCetsUnknown()->getIdentifier())
  );
}
 
void
PacemakerBidomainProblem::run()
{
  if (ceps::isProfiling())
    getProfiler()->start("whole", "solving the bidomain problem");
 
  initializeEquation();
 
  // Solve the static problem
  if (m_parameters->isActiveOption("compute electric field"))
  {
   auto spb = ceps::getNew<FluxAnodeCathodeProblem>(m_geom, m_parameters);
    auto dico = spb->getFunctionDictionary();
    dico->deleteTensor("physCoeffSigma");
    dico->add("sigma_i", m_sigmaIOptions);
    dico->add("sigma_e", m_sigmaEOptions);
    dico->add("physCoeffSigma", "OPERATOR sigma_i + sigma_e");
    spb->setOutputFileBase(spb->getOutputFileBase() + "_elec_field");
    spb->initializeEquation();
 
    FluxAnodeCathodeSolver ssolver(spb);
    ssolver.solve();
    ceps::destroyObject(spb);
  }
 
  PacemakerBidomainSolver solver(this);
 
  solver.solve();
 
  if (hasAnalyticSolution())
    m_errors = solver.getErrors();
 
  if (ceps::isProfiling())
    getProfiler()->stop("whole");
}
 
CepsSet<CepsAttribute>
PacemakerBidomainProblem::getAnodeAttributes() const
{
  return getAnode()->getAttributes();
}
 
CepsSet<CepsAttribute>
PacemakerBidomainProblem::getCathodeAttributes() const
{
  return getCathode()->getAttributes();
}
 
Unknown*
PacemakerBidomainProblem::getUAnodeUnknown() const
{
  return getUnknown(UnknownsName::UAnode);
}
 
Unknown*
PacemakerBidomainProblem::getUCathodeUnknown() const
{
  return getUnknown(UnknownsName::UCathode);
}
 
Unknown*
PacemakerBidomainProblem::getVCetsUnknown() const
{
  return getUnknown(UnknownsName::VCets);
}
 
Unknown*
PacemakerBidomainProblem::getVCtsUnknown() const
{
  return getUnknown(UnknownsName::VCts);
}
 
void
PacemakerBidomainProblem::setupWithParameters(InputParameters* params)
{
  // enforce no anodal and cathodal stimulation for surrounded bidomain problem
  m_biElecStimOpts = "";
 
  return;
}