CardiacProblem.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/CardiacProblem.hpp"
#include "cardiac/solver/CardiacSolver.hpp"
 
// -----------------------------------------------------------------------------------------------
// PUBLIC
 
CardiacProblem::
CardiacProblem(Geometry* g, InputParameters* params): 
  HeatProblem         (g,params),
  m_CmOptions         (""),
  m_AmOptions         ("CONSTANT 1000"),
  m_ionOptions        (""),
  m_fibersOptions     (""),
  m_sigmaIOptions     ("CONSTANT 0.4 0.2 0.2"),
  m_sigmaEOptions     ("CONSTANT 0.4 0.2 0.2"),
  m_volFracOptions    (""),
  m_volFracItpFileI   (""),
  m_volFracItpFileE   (""),
  m_stimOptions       (""),
  m_bcStimOptions     (""),
  m_tissueAttributes  ({-1}),
  m_Cm                (nullptr),
  m_Am                (nullptr),
  m_ionicModels       ({}),
  m_sigmaI            (nullptr),
  m_sigmaE            (nullptr),
  m_volFrac           (nullptr)
{
  if (ceps::isValidPtr(params))
    CardiacProblem::setupWithParameters(params);
}
 
CardiacProblem::
~CardiacProblem()
{
  ceps::destroyObject(m_Cm);
  ceps::destroyObject(m_Am);
  for (auto ionM: m_ionicModels)
    ceps::destroyObject(ionM);
  ceps::destroyObject(m_sigmaI);
  ceps::destroyObject(m_sigmaE);
  ceps::destroyObject(m_volFrac);
}
 
void
CardiacProblem::getInitialCondition(DHVecPtr vm) const
{
  vm->zero();
  for (auto model: m_ionicModels)
    model->getLinkedSolver()->fillPdeVm(vm);
}
 
CepsVector<AbstractIonicModel*>
CardiacProblem::getIonicModels() const
{
  return m_ionicModels;
}
 
ScalarField<DegreeOfFreedom>*
CardiacProblem::getAm() const
{
  return m_Am;
}
 
ScalarField<DegreeOfFreedom>*
CardiacProblem::getCm() const
{
  return m_Cm;
}
 
TensorField<GeomCell>*
CardiacProblem::getSigmaI() const
{
  return m_sigmaI;
}
 
TensorField<GeomCell>*
CardiacProblem::getSigmaE() const
{
  return m_sigmaE;
}
 
const CepsSet<CepsAttribute>&
CardiacProblem::getTissueAttributes() const
{
  return m_tissueAttributes;
}
 
 
// -------------------------------------------------------------------------------------------------
// Protected
 
void
CardiacProblem::initializeEquation()
{
  this->defineUnknowns();
  this->createSpatialDiscretization();
 
  CEPS_SEPARATOR;
  CEPS_SAYS("Setting up elements of the cardiac problem");
  initializeIonicModels();
  initializeCm();
  initializeAm();
  this->initializeConductivities();
 
  this->defineSourceTerms();
  this->defineBoundaryConditions();
  this->defineAnalyticSolution();
}
 
void 
CardiacProblem::setupWithParameters(InputParameters* params)
{
 
  CepsSet<CepsInt> attrs = ceps::toInts(params->getString("tissue attributes",""));
  if (not attrs.empty())
  {
    m_tissueAttributes.clear();
    m_tissueAttributes.insert(attrs.begin(),attrs.end());
  }
 
  m_ionOptions = params->getString("ionic model");
 
  if (params->hasKey("membrane capacitance"))
    m_CmOptions = params->getString("membrane capacitance",m_CmOptions);
  else
    m_CmOptions = params->getString("cm",m_CmOptions);
 
  if (params->hasKey("membrane surface"))
    m_AmOptions = params->getString("membrane surface",m_AmOptions);
  else
    m_AmOptions = params->getString("am",m_AmOptions);
 
  m_fibersOptions = params->getString("fibers",m_fibersOptions);
 
  if (params->hasKey("sigma i"))
    m_sigmaIOptions = params->getString("sigma i",m_sigmaIOptions);
  else
    m_sigmaIOptions = params->getString("intracellular conductivity",m_sigmaIOptions);  
  if (params->hasKey("sigma e"))
    m_sigmaEOptions = params->getString("sigma e",m_sigmaEOptions);
  else
    m_sigmaEOptions = params->getString("extracellular conductivity",m_sigmaEOptions);  
 
  m_volFracOptions  = params->getString("volume fraction",m_volFracOptions);
  m_volFracItpFileI = params->getString("volume fraction map intracellular",m_volFracItpFileI);
  m_volFracItpFileE = params->getString("volume fraction map extracellular",m_volFracItpFileE);
 
  m_stimOptions   = params->getString("stimulation"         ,m_stimOptions  );
  m_bcStimOptions = params->getString("boundary stimulation",m_bcStimOptions);
}
 
void
CardiacProblem::initializeCm()
{
 
  if (not m_CmOptions.empty())
    m_functions->add("dicoCm",m_CmOptions);
  else
  // Build Cm vector from ionic models in case it is needed by the PDE.
  {
    CepsMap<CepsAttribute,CepsReal> cmMap;
    for (auto model: m_ionicModels)
      for (auto attr: model->getAttributes())
        cmMap[attr] = model->getPaperCm(true);
    // Add a -1 entry in case it is not defined everywhere
    // it is a default value : eval starts with meaningful attributes
    cmMap[-1] = 1.;
    m_functions->addEntry("dicoCm",newCstPiecewiseSAFunc(cmMap));
  }
 
  auto dofs = m_discr->getDistributedDofs();
  m_Cm = ceps::getNew<ScalarField<DegreeOfFreedom>>(m_functions->getScalar("dicoCm"),&(dofs->getOwned()));
 
  // Link Cm built from options to ionic models
  if (not m_CmOptions.empty())
    for (auto model: m_ionicModels)
      model->setCm(m_Cm);
 
  return;
}
 
void
CardiacProblem::initializeAm()
{
  m_functions->add("dicoAm",m_AmOptions);
  auto dofs = m_discr->getDistributedDofs();
  m_Am = ceps::getNew<ScalarField<DegreeOfFreedom>>(m_functions->getScalar("dicoAm"),&(dofs->getOwned()));
  return;
}
 
void
CardiacProblem::initializeIonicModels()
{
 
  m_ionicModels.clear();
 
  CEPS_ABORT_IF(m_ionOptions.empty(),
    "There is no ionic model provided as input !"
  );
  // This map keeps track of all ionic models region assignments for all unknowns
  CepsMap<CepsUnknownIndex,CepsSet<CepsAttribute>> prevAttrs;
  auto tAttrs = this->getTissueAttributes();
 
  for (CepsString ionOpt : ceps::split(m_ionOptions,","))
    m_ionicModels.push_back(getIonicModelFromString(ionOpt,m_parameters,this,prevAttrs,tAttrs));
   
  return;
 
}
 
void
CardiacProblem::initializeConductivities()
{
  this->initializeConductivities(m_sigmaIOptions,m_sigmaEOptions,m_sigmaI,m_sigmaE,m_fibersOptions);
}
 
void
CardiacProblem::initializeConductivities(
  CepsString              sigmaIOpts, 
  CepsString              sigmaEOpts,
  TensorField<GeomCell>*& sigmaI, 
  TensorField<GeomCell>*& sigmaE,
  CepsString              fiberOptions,
  CepsString              suffix
)
{
  CepsBool doVolFrac = not m_volFracOptions.empty();
 
  // ==========================================================
  // Volume fraction : a reference tensor is interpolated from file 
  if (doVolFrac)
  {
    auto words = ceps::split(m_volFracOptions);
 
    CepsString errMessage = 
      CepsString("Incorrect options for volume fraction: please provide either :\n") +
      "     - a file name and the name of the scalar data array, defined on cells.\n" +
      "     - a string defining a physical coefficient (cf documentation)\n";
 
    CEPS_ABORT_IF(words.size()<1,errMessage);
 
    if (ceps::fileExists(words[0]))
    {
      CEPS_ABORT_IF(words.size()<2,errMessage);
      m_functions->add("volFrac","FILE " + words[0] + " SCALAR " + words[1] + " CELL",getGeometry());
    }
    else
    {
      m_functions->add("volFrac",m_volFracOptions);
    }
 
    // Create the scalar field with vol frac
 
    // Create interpolation table
    CepsMap<CepsReal,CepsMathTensor> tableI,tableE;
    std::ifstream ifileI(m_volFracItpFileI);
    CEPS_ABORT_IF(not ifileI.good(),
      "Unable to open file with volume fraction interpolation map : " << m_volFracItpFileI
    );
 
    CepsUInt n = ceps::readInt(ifileI,"Unable to read number of lines in volume fraction table file");
    for (CepsUInt i=0;i<n;i++)
    {
      CepsReal       x = ceps::readReal  (ifileI,"Cannot read coefficient in volume fraction table");
      CepsMathTensor t = ceps::readTensor(ifileI,"Cannot read tensor in volume fraction table");
      tableI[x] = t;
    }
    std::ifstream ifileE(m_volFracItpFileE);
    CEPS_ABORT_IF(not ifileE.good(),
      "Unable to open file with volume fraction interpolation map : " << m_volFracItpFileE
    );
 
    n = ceps::readInt(ifileE,"Unable to read number of lines in volume fraction table file");
    for (CepsUInt i=0;i<n;i++)
    {
      CepsReal       x = ceps::readReal  (ifileE,"Cannot read coefficient in volume fraction table");
      CepsMathTensor t = ceps::readTensor(ifileE,"Cannot read tensor in volume fraction table");
      tableE[x] = t;
    }
 
    ScalarSAFunc* funcV = ceps::runtimeCast<ScalarSAFunc*>(m_functions->getScalar("volFrac"));
 
    m_functions->addEntry("sigmaIVF"+suffix,newCoeffInterpolatorSAFunc(tableI,funcV));
    m_functions->addEntry("sigmaEVF"+suffix,newCoeffInterpolatorSAFunc(tableE,funcV));
 
    // We modify conductivity only in the tissue !
    m_functions->add("sigmaIOrg"+suffix,sigmaIOpts);
    m_functions->add("sigmaEOrg"+suffix,sigmaEOpts);
 
    for (CepsString tag: {"I","E"})
    {
      CepsString tmpopts = "FPIECEWISE "; 
      for (CepsAttribute attr : this->getTissueAttributes())
        tmpopts += ceps::toString(attr) + " sigma"+tag+"VF"+suffix+", ";
      tmpopts += "-1 sigma"+tag+"Org"+suffix;
      m_functions->add("sigma"+tag+"Ref"+suffix,tmpopts);
    }
 
  }
  // ==========================================================
  // No volume fraction : a reference diagonal tensor is read from options
  else
  {
    m_functions->add("sigmaIRef"+suffix,sigmaIOpts);
    m_functions->add("sigmaERef"+suffix,sigmaEOpts);
  }
 
  createFibersAndConductivities(fiberOptions,sigmaI,sigmaE,suffix);
 
}
 
void
CardiacProblem::initializeStimulations()
{
 
  CepsUInt i=0;
  for (CepsString opts: ceps::split(m_stimOptions,","))
  {
    // Add function to dictionary
    CepsString label = "stim"+ceps::toString(i);
    m_functions->addFunction("dico"+label,opts);
 
    CepsString fullOpts = "dico" + label + " " + opts + " TYPE IAPP";
    if (this->getTMVUnknowns().size()==1 and opts.find("UNKNOWN") == opts.npos)
      fullOpts += " UNKNOWN 0";
 
    m_sourceTerms->add(label,fullOpts,true);
 
    // Link stimulation to corresponding ionic model
    ScalarSourceTerm* stim = m_sourceTerms->getSourceTerm(label); // we just created it
    Unknown*          u    = this->getUnknown(stim->getUnknownId());
    for (AbstractIonicModel* model: m_ionicModels)
      if (model->getUnknown() == u)
        model->addStimulation(stim);
 
    i++;
  }
 
}
 
void
CardiacProblem::defineBoundaryConditions()
{
  CepsUInt i=0;
  for (CepsString opts: ceps::split(m_bcStimOptions,","))
  {
    // Add function to dictionary
    CepsString label = "bcstim"+ceps::toString(i);
    m_functions->addFunction("dico"+label,opts);
 
    CepsString fullOpts = "dico" + label + " " + opts;
    
    if (this->getTMVUnknowns().size()==1 and opts.find("UNKNOWN") == opts.npos)
      fullOpts += " UNKNOWN 0";
 
    if (opts.find("DIRICHLET") != opts.npos)
      m_boundaryConditions->add(label,CepsBoundaryConditionFlag::Dirichlet,fullOpts);
    else if (opts.find("NEUMANN") != opts.npos)
      m_boundaryConditions->add(label,CepsBoundaryConditionFlag::Neumann,fullOpts);
    else
    {
      CEPS_ABORT(
        "When specifying a boundary stimulation, the DIRICHLET or "
        "NEUMANN keyword must be given in the options"
      );
    }
 
    i++;
  }
 
}
 
void
CardiacProblem::defineSourceTerms()
{
  initializeStimulations();
}
 
void
CardiacProblem::createFibersAndConductivities(
  CepsString              fiberOptions,
  TensorField<GeomCell>*& sigmaI,
  TensorField<GeomCell>*& sigmaE,
  CepsString              suffix
)
{
  CepsBool doVolFrac = not m_volFracOptions.empty();
  
  if (fiberOptions.empty())
  {
    m_functions->add("fibersL"+suffix,"CONSTANT 1. 0. 0.");
    m_functions->add("fibersT"+suffix,"CONSTANT 0. 1. 0.");
    m_functions->add("fibersN"+suffix,"CONSTANT 0. 0. 1.");
  }
  else
  {
    auto words = ceps::split(fiberOptions);
 
    if (ceps::fileExists(words[0]))
    {
      CepsString keywordL = words.size() < 2 ? "fibers"            : words[1];
      CepsString keywordT = words.size() < 3 ? "fibers_transverse" : words[2];
      CepsString keywordN = words.size() < 4 ? "fibers_normal"     : words[3];
      m_functions->addFromFile("fibersL"+suffix,words[0]+" VECTOR "+keywordL+" CELL",getGeometry());
      m_functions->addFromFile("fibersT"+suffix,words[0]+" VECTOR "+keywordT+" CELL",getGeometry(),true);
      m_functions->addFromFile("fibersN"+suffix,words[0]+" VECTOR "+keywordN+" CELL",getGeometry(),true);
    }
    else
    {
      auto fopts = ceps::split(fiberOptions,",");
      m_functions->add("fibersL"+suffix,fopts[0],m_geom);
      if (fopts.size()>1)
        m_functions->add("fibersT"+suffix,fopts[1],m_geom);
      if (fopts.size()>2)
        m_functions->add("fibersN"+suffix,fopts[2],m_geom);
    }
 
  }
 
  VectorSAFunc* fL = m_functions->getVector("fibersL"+suffix);
  VectorSAFunc* fT = m_functions->getVector("fibersT"+suffix,nullptr);
  VectorSAFunc* fN = m_functions->getVector("fibersN"+suffix,nullptr);
 
  CEPS_ABORT_IF(doVolFrac and ceps::isNullPtr(fT),
    "Using volume fraction requires that at least two components of the fibers orientation are provided"
  );
 
  if (m_functions->hasTensor("sigmaIRef"))
  {
    m_functions->addEntry("sigmaI"+suffix,new FiberAligner(fL,fT,fN,m_functions->getTensor("sigmaIRef")));
    m_functions->addEntry("sigmaE"+suffix,new FiberAligner(fL,fT,fN,m_functions->getTensor("sigmaERef")));
  }
  else
  {
    m_functions->addEntry("sigmaI"+suffix,new FiberAligner(fL,fT,fN,nullptr,m_functions->getVector("sigmaIRef")));
    m_functions->addEntry("sigmaE"+suffix,new FiberAligner(fL,fT,fN,nullptr,m_functions->getVector("sigmaERef")));
  }
 
  // ==========================================================
  // Create fields
 
  // We do not bufferize data, unless there is volume fraction, which makes conductivity
  // different at each point
  CepsVector<GeomCell*> emp;
  sigmaI = ceps::getNew<TensorField<GeomCell>>(m_functions->getTensor("sigmaI"+suffix),&emp,doVolFrac);
  sigmaE = ceps::getNew<TensorField<GeomCell>>(m_functions->getTensor("sigmaE"+suffix),&emp,doVolFrac);
 
  // In optimized mode adding the empty vector domain does not seem to work
  // size is undefined, very large
  sigmaI->clearDomains();
  sigmaE->clearDomains();
 
  for(CepsUInt dim=1;dim<=3;dim++)
    if (m_geom->hasMeshOfDim(dim))
    {
      sigmaI->addDomain(&(m_geom->getMeshOfDim(dim)->getCells()));
      sigmaE->addDomain(&(m_geom->getMeshOfDim(dim)->getCells()));
    }
 
}
 
CardiacProblem::FiberAligner::FiberAligner(
  VectorSAFunc* fL, 
  VectorSAFunc* fT, 
  VectorSAFunc* fN,
  TensorSAFunc* sT,
  VectorSAFunc* sV
):
  m_fL(fL),m_fT(fT),m_fN(fN),m_sT(sT),m_sV(sV)
{
}
 
CepsMathTensor
CardiacProblem::FiberAligner::eval(CepsStandardArgs args) 
{
 
  if (ceps::isNullPtr(m_fT)) // one direction
  {
    CepsReal       g0    = ceps::isNullPtr(m_sT) ? m_sV->eval(args)(0) : m_sT->eval(args)(0,0);
    CepsReal       g1    = ceps::isNullPtr(m_sT) ? m_sV->eval(args)(1) : m_sT->eval(args)(1,1);
    CepsMathVertex f     = m_fL->eval(args);
    CepsMathTensor outer = f*f.transpose();
    return g1*CepsMathTensor::Identity() + (g0-g1)*outer;
  }
 
  CepsMathTensor f;
  CepsMathTensor s;
  if (ceps::isNullPtr(m_sT))
    s = m_sV->eval(args).asDiagonal();
  else
    s = m_sT->eval(args);
 
  f.col(0) = m_fL->eval(args);
  f.col(1) = m_fT->eval(args);
  if (ceps::isNullPtr(m_fN)) // two directions
    f.col(2) = f.col(0).cross(f.col(1));
  else
    f.col(2) = m_fN->eval(args);
 
  return f*s*f.transpose();
 
}
 
CepsEnum
CardiacProblem::FiberAligner::getFlags() const
{
  CepsEnum res1 = m_fL->getFlags();
  if (ceps::isValidPtr(m_fT))
    res1 = res1 | m_fT->getFlags();
  if (ceps::isValidPtr(m_fN))
    res1 = res1 | m_fN->getFlags();
 
  CepsEnum res2 = ceps::isValidPtr(m_sT) ? m_sT->getFlags() : m_sV->getFlags();
 
  return res1 | res2;
 
}