PacemakerPoissonProblem.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/PacemakerPoissonProblem.hpp"
 
#include "cardiac/pacemaker/solver/PacemakerPoissonSolver.hpp"
#include "pde/problem/static/FluxAnodeCathodeProblem.hpp"
#include "pde/solver/static/FluxAnodeCathodeSolver.hpp"
 
PacemakerPoissonProblem::PowerSupply::PowerSupply(PacemakerPoissonProblem* problem) :
    AbstractPacemakerProblem::PowerSupply(problem)
{}
 
CepsReal
PacemakerPoissonProblem::PowerSupply::eval(CepsStandardArgs args)
{
  PacemakerPoissonProblem* pb = ceps::runtimeCast<PacemakerPoissonProblem*>(m_problem);
 
  if (args.unknownId == pb->getVCetsUnknown()->getIdentifier())
    return AbstractPacemakerProblem::PowerSupply::eval(args);
  return 0.0;
}
 
PacemakerPoissonProblem::PacemakerPoissonProblem(Geometry* g, InputParameters* params) :
    HeatProblem(g, params),
    AbstractPacemakerProblem(g, params)
{
  setProblemName("Poisson connected to a Pacemaker");
 
  if (ceps::isValidPtr(params))
    PacemakerPoissonProblem::setupWithParameters(params);
}
 
void
PacemakerPoissonProblem::getInitialCondition(DHVecPtr v0) const
{
  v0->zero();
  HeatProblem::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
PacemakerPoissonProblem::defineUnknowns()
{
  addUnknown(UnknownsName::u, {CepsUniversal});  // ue is defined everywhere
 
  // the 0D unknowns "potential"
  addZeroDUnknown(UnknownsName::UAnode);
  addZeroDUnknown(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);
  addUnknownInteraction(UnknownsName::UAnode, UnknownsName::UCathode);
 
  addUnknownInteraction(UnknownsName::UCathode, UnknownsName::u, getCathode()->getAttributes());
  addUnknownInteraction(UnknownsName::UAnode, UnknownsName::u, getCathode()->getAttributes());
}
 
void
PacemakerPoissonProblem::defineBoundaryConditions()
{}
 
void
PacemakerPoissonProblem::defineSourceTerms()
{
  HeatProblem::defineSourceTerms();
 
  m_functions->addEntry("capaChargeFunc", ceps::getNew<PowerSupply>(this));
  m_sourceTerms->add(
    "__power_supply capaChargeFunc UNKNOWN "  //
    + ceps::toString(getVCetsUnknown()->getIdentifier())
  );
}
 
void
PacemakerPoissonProblem::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);
    spb->getFunctionDictionary()->deleteTensor("physCoeffSigma");
    spb->getFunctionDictionary()->addConstant("physCoeffSigma", m_k);
    spb->setOutputFileBase(spb->getOutputFileBase() + "_elec_field");
    spb->initializeEquation();
 
    FluxAnodeCathodeSolver ssolver(spb);
    ssolver.solve();
    ceps::destroyObject(spb);
  }
 
  PacemakerPoissonSolver solver(this);
 
  solver.solve();
 
  if (hasAnalyticSolution())
    m_errors = solver.getErrors();
 
  if (ceps::isProfiling())
    getProfiler()->stop("whole");
}
 
CepsSet<CepsAttribute>
PacemakerPoissonProblem::getAnodeAttributes() const
{
  return getAnode()->getAttributes();
}
 
CepsSet<CepsAttribute>
PacemakerPoissonProblem::getCathodeAttributes() const
{
  return getCathode()->getAttributes();
}
 
Unknown*
PacemakerPoissonProblem::getUFieldUnknown() const
{
  return getUnknown(UnknownsName::u);
}
 
Unknown*
PacemakerPoissonProblem::getUAnodeUnknown() const
{
  return getUnknown(UnknownsName::UAnode);
}
 
Unknown*
PacemakerPoissonProblem::getUCathodeUnknown() const
{
  return getUnknown(UnknownsName::UCathode);
}
 
Unknown*
PacemakerPoissonProblem::getVCetsUnknown() const
{
  return getUnknown(UnknownsName::VCets);
}
 
Unknown*
PacemakerPoissonProblem::getVCtsUnknown() const
{
  return getUnknown(UnknownsName::VCts);
}
 
void
PacemakerPoissonProblem::setupWithParameters(InputParameters* params)
{
  return;
}