BilayerMonodomainProblem.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 "BilayerMonodomainProblem.hpp"
#include "cardiac/solver/BilayerMonodomainSolver.hpp"
 
BilayerMonodomainProblem::BilayerMonodomainProblem(Geometry* g, InputParameters* params): 
  CardiacProblem(g,params),
  m_couplingStr    (200.),
  m_couplingAttrs  ({}),
  m_couplingAttrNeg(false),
  m_sigmaIOptions2 ("CONSTANT 1. 1. 1."),
  m_sigmaEOptions2 ("CONSTANT 1. 1. 1."),
  m_fibersOptions2 (""),
  m_sigmaI2        (nullptr),
  m_sigmaE2        (nullptr)
{
  setProblemName("Bilayer Monodomain");
 
  if(ceps::isValidPtr(params))
    BilayerMonodomainProblem::setupWithParameters(params);
 
  CEPS_SAYS( 
    "Bilayer model: elements " << (m_couplingAttrNeg ? "without":"with") << 
    " attribute " << m_couplingAttrs << " will be coupled,\n" << 
    "adding the equality " << m_couplingStr << "(u1-u2) = 0 to the original equation."
  );
}
 
void
BilayerMonodomainProblem::defineUnknowns()
{
  addUnknown(UnknownsName::vm1, m_tissueAttributes);
  addUnknown(UnknownsName::vm2, m_tissueAttributes);
}
 
CepsVector<Unknown*>
BilayerMonodomainProblem::getTMVUnknowns() const
{
  return {getTMVUnknown1(), getTMVUnknown2()};
}
 
CepsVector<Unknown*>
BilayerMonodomainProblem::getCardiacUnknowns() const
{
  return getTMVUnknowns();
}
 
Unknown*
BilayerMonodomainProblem::getTMVUnknown1() const
{
  return getUnknown(UnknownsName::vm1);
}
 
Unknown*
BilayerMonodomainProblem::getTMVUnknown2() const
{
  return getUnknown(UnknownsName::vm2);
}
 
TensorField<GeomCell>*
BilayerMonodomainProblem::getSigmaI2() const
{
  return m_sigmaI2;
}
 
TensorField<GeomCell>*
BilayerMonodomainProblem::getSigmaE2() const
{
  return m_sigmaE2;
}
 
void
BilayerMonodomainProblem::run()
{
  if (ceps::isProfiling()) 
    getProfiler()->start("whole","performing whole simulation (BilayerMonodomainProblem)");
 
  initializeEquation();
 
  BilayerMonodomainSolver solver(this);
  
  solver.solve(); 
 
  if (hasAnalyticSolution())
    m_errors = solver.getErrors();
 
  if (ceps::isProfiling())
    getProfiler()->stop("whole");  
 
}
 
CepsReal
BilayerMonodomainProblem::getCouplingStrength() const
{
  return m_couplingStr;
}
 
CepsSet<CepsAttribute>
BilayerMonodomainProblem::getCouplingAttributes() const
{
  return m_couplingAttrs;
}
 
CepsBool
BilayerMonodomainProblem::getCouplingAttrNeg() const
{
  return m_couplingAttrNeg;
}
 
void
BilayerMonodomainProblem::initializeConductivities()
{ 
  // Creates reference sigmaI and sigmaE, 
  // then reorient sigmaI and sigmaE for first layer
  CardiacProblem::initializeConductivities(
    m_sigmaIOptions,
    m_sigmaEOptions,
    m_sigmaI,
    m_sigmaE,
    m_fibersOptions
  );
  // Second layer
  CardiacProblem::initializeConductivities(
    m_sigmaIOptions2,
    m_sigmaEOptions2,
    m_sigmaI2,
    m_sigmaE2,
    m_fibersOptions2,
    "2"
  );
 
}
 
void
BilayerMonodomainProblem::setupWithParameters(InputParameters* params)
{
 
  // Conductivities and fibers
  for (CepsVector<CepsString> v: CepsVector<CepsVector<CepsString>>(
    {{m_sigmaIOptions ,m_sigmaEOptions ,"1"},
     {m_sigmaIOptions2,m_sigmaEOptions2,"2"}}
  ))
  {
    if (params->hasKey("sigma i layer "+v[2]))
      m_sigmaIOptions = params->getString("sigma i layer "+v[2],v[0]);
    else
      m_sigmaIOptions = params->getString("intracellular conductivity layer "+v[2],v[0]);  
    if (params->hasKey("sigma e layer "+v[2]))
      m_sigmaEOptions = params->getString("sigma e layer "+v[2],v[1]);
    else
      m_sigmaEOptions = params->getString("extracellular conductivity layer "+v[2],v[1]);  
  }
 
  m_fibersOptions  = params->getString("fibers layer 1",m_fibersOptions2);
  m_fibersOptions2 = params->getString("fibers layer 2",m_fibersOptions2);
 
  // Coupling
  m_couplingStr = params->getReal("coupling strength",m_couplingStr);
 
  CepsVector<CepsString> ws = ceps::split(params->getString("coupling attribute"));
  m_couplingAttrNeg = ceps::toUpper(ws[0])=="NOT";
  if (m_couplingAttrNeg)
  {
    CEPS_ABORT_IF(ws.size()<2,
      "BilayerParameters: could not read attribute of non-coupled elements"
    );
    for(CepsUInt i=1;i<ws.size();i++) 
      m_couplingAttrs.insert(ceps::toInt(ws[i]));
  }
  else
  {
    for(CepsUInt i=0;i<ws.size();i++) 
      m_couplingAttrs.insert(ceps::toInt(ws[i]));
    CepsBool ok = m_tissueAttributes.empty() or m_tissueAttributes.contains(-1)
                or m_couplingAttrs.contains(-1);
 
    for (CepsAttribute attr:m_couplingAttrs)
      ok = ok or m_tissueAttributes.contains(attr);
 
    CEPS_ABORT_IF(not ok, 
      "Any of the given coupling attributes can not be found in the tissue attributes"
    );
  }
 
}