ExtendedBidomainProblem.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/problem/ExtendedBidomainProblem.hpp"
 
#include "cardiac/solver/ExtendedBidomainSolver.hpp"
#include "pde/problem/static/FluxAnodeCathodeProblem.hpp"
#include "pde/solver/static/FluxAnodeCathodeSolver.hpp"
 
ExtendedBidomainProblem::ExtendedBidomainProblem(Geometry* g, InputParameters* params) :
    BidomainProblem(g, params),
    m_extMediaAttrsOpts(""),
    m_extMediaAttrs({}),
    m_biElecStimOpts(""),
    m_robinCoeffs(""),
    m_anodeAttrsOpts(""),
    m_anodeAttrs({}),
    m_cathodeAttrsOpts(""),
    m_cathodeAttrs({})
{
  setProblemName("Extended Bidomain");
 
  if (ceps::isValidPtr(params))
    ExtendedBidomainProblem::setupWithParameters(params);
 
  // Set electrode attributes
  m_anodeAttrs   = ceps::toInts(m_anodeAttrsOpts);
  m_cathodeAttrs = ceps::toInts(m_cathodeAttrsOpts);
 
  m_extMediaAttrs = ceps::toInts(m_extMediaAttrsOpts);
  if (m_extMediaAttrs.empty())
    m_extMediaAttrs = {-1};
 
  m_extMediaAttrs.insert(m_anodeAttrs.begin(), m_anodeAttrs.end());
  m_extMediaAttrs.insert(m_cathodeAttrs.begin(), m_cathodeAttrs.end());
 
  // CEPS_ABORT_IF(
  //   m_extMediaAttrs.empty(),
  //   "missing blood/torso/extracardiac attributes to run an extended bidomain model"
  // );
}
 
void
ExtendedBidomainProblem::defineUnknowns()
{
  addUnknown(UnknownsName::vm, getTissueAttributes());
  addUnknown(UnknownsName::ue, getWholeDomainAttributes()); // ue is defined everywhere
 
  addUnknownInteraction(UnknownsName::vm, UnknownsName::ue, getTissueAttributes());
 
  if (requireNullMean())
  {
    // and a lagrangian unknown
    addZeroDUnknown(UnknownsName::lag);
    addUnknownInteraction(UnknownsName::lag, UnknownsName::ue);
  }
}
 
 
void
ExtendedBidomainProblem::run()
{
  if (ceps::isProfiling())
    getProfiler()->start("whole", "solving the bidomain problem");
 
  initializeEquation();
 
  // Solve the static problem
  if (hasAnodalCathodalStimulation() and m_parameters->isActiveOption("compute electric field"))
  {
    FluxAnodeCathodeProblem spb(m_geom, m_parameters);
    spb.getFunctionDictionary()->deleteTensor("physCoeffSigma");
    spb.getFunctionDictionary()->add("sigma_i", m_sigmaIOptions);
    spb.getFunctionDictionary()->add("sigma_e", m_sigmaEOptions);
    spb.getFunctionDictionary()->add("physCoeffSigma", "OPERATOR sigma_i + sigma_e");
    spb.setOutputFileBase(spb.getOutputFileBase() + "_elec_field");
    spb.initializeEquation();
 
    FluxAnodeCathodeSolver ssolver(&spb);
    ssolver.solve();
  }
 
  ExtendedBidomainSolver solver(this);
 
  solver.solve();
 
  if (hasAnalyticSolution())
    m_errors = solver.getErrors();
 
  if (ceps::isProfiling())
    getProfiler()->stop("whole");
}
 
const CepsSet<CepsAttribute>&
ExtendedBidomainProblem::getExtracardiacAttributes() const
{
  return m_extMediaAttrs;
}
 
CepsSet<CepsAttribute>
ExtendedBidomainProblem::getWholeDomainAttributes() const
{
  auto tissue = getTissueAttributes();
  auto extra  = getExtracardiacAttributes();
 
  CepsSet<CepsAttribute> res = {};
  res.insert(tissue.begin(),tissue.end());
  res.insert(extra .begin(),extra .end());
  return res;
}
 
CepsSet<CepsAttribute>
ExtendedBidomainProblem::getAnodeAttributes() const
{
  return m_anodeAttrs;
}
 
CepsSet<CepsAttribute>
ExtendedBidomainProblem::getCathodeAttributes() const
{
  return m_cathodeAttrs;
}
 
CepsBool
ExtendedBidomainProblem::hasAnodalCathodalStimulation() const
{
  CepsBool res = true;
  res          = res and not m_biElecStimOpts.empty();
  res          = res and m_anodeAttrs  .size() > 0;
  res          = res and m_cathodeAttrs.size() > 0;
  return res;
}
 
void
ExtendedBidomainProblem::defineBoundaryConditions()
{
  if (hasAnodalCathodalStimulation())
  {
    
    m_functions->add("bielec_stim", CepsString("FUNCTION TIME ") + m_biElecStimOpts);
 
    m_functions->add("anode_stim", "OPERATOR 1.0 * bielec_stim");
    m_functions->add("cathode_stim", "OPERATOR -1.0 * bielec_stim");
 
    m_boundaryConditions->add(
      "anode NEUMANN anode_stim "
      "UNKNOWN 1 ATTRIBUTES " +
      m_anodeAttrsOpts
    );
    m_boundaryConditions->add(
      "cathode NEUMANN cathode_stim "
      "UNKNOWN 1 ATTRIBUTES " +
      m_cathodeAttrsOpts
    );
     
    for (auto definition : ceps::split(m_robinCoeffs, ","))
    {
      auto options = ceps::split(definition);
      CEPS_ABORT_IF(options.size() != 2, "robin coefficient entry requires a scalar and an attribute");
      m_boundaryConditions->add("boundary ROBIN 0.0 UNKNOWN 1 ATTRIBUTES " + options[1] + " ROBINCOEFF " + options[0]);
    }
 
   /* 
    m_boundaryConditions->add(
      "anode DIRICHLET bielec_stim "
      "UNKNOWN 1 ATTRIBUTES " +
      m_anodeAttrsOpts
    );
    m_boundaryConditions->add(
      "cathode DIRICHLET 0. "
      "UNKNOWN 1 ATTRIBUTES " +
      m_cathodeAttrsOpts
    );
    */
 
 
 
  }
}
 
CepsBool
ExtendedBidomainProblem::requireNullMean() const
{
  return m_robinCoeffs.empty();
} 
 
void
ExtendedBidomainProblem::setupWithParameters(InputParameters* params)
{
  // blood medium
  m_extMediaAttrsOpts += params->getString("blood attributes", "");
  // torso attributes, large domain usually
  m_extMediaAttrsOpts += " " + params->getString("torso attributes", "");
  // generic extracardiac attributes
  m_extMediaAttrsOpts += " " + params->getString("extracardiac attributes", "");
 
  // anodal and cathodal stimulation
  m_biElecStimOpts = params->getString("electrodes stimulation", m_biElecStimOpts);
 
  m_anodeAttrsOpts   = params->getString("anode attributes", m_anodeAttrsOpts);
  m_cathodeAttrsOpts = params->getString("cathode attributes", m_cathodeAttrsOpts);
 
  m_robinCoeffs = params->getString("robin coefficients", "");
 
  return;
}