CLMonodomainProblem.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/CLMonodomainProblem.hpp"
#include "cardiac/solver/MonodomainSolver.hpp"
#include "pde/problem/static/FluxAnodeCathodeProblem.hpp"
#include "pde/solver/static/FluxAnodeCathodeSolver.hpp"
 
 
CLMonodomainProblem::Ufunc::Ufunc(DHVecPtr* sol, CepsFunctionFlag loc) :
  SolVecSAFunc(sol,loc)
{
}
 
CepsReal
CLMonodomainProblem::Ufunc::eval(CepsStandardArgs args)
{
  setArgsFromLocation(args);
  return args.u;
}
 
CLMonodomainProblem::CLMonodomainProblem(Geometry* g, InputParameters* params)
  : MonodomainProblem(g,params)
{
  setProblemName("Current Lifted Monodomain");
}
 
void
CLMonodomainProblem::run()
{
  if (ceps::isProfiling())
    getProfiler()->start("whole","solving the whole PDE (FluxAnodeCathodeProblem + CLMonodomainProblem)");
 
  initializeEquation();
 
  // Solve the static problem
 
  FluxAnodeCathodeProblem spb(m_geom, m_parameters);
  // if (spb.getFunctionDictionary()->hasTensor("physCoeffSigma"))
  //   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() + "_static");
  spb.initializeEquation();
 
  FluxAnodeCathodeSolver ssolver(&spb);
  ssolver.solve(); 
  
  // not realy clean, but it's working
  auto solCL = ssolver.getSolution();
  getFunctionDictionary()->addEntry("ubar", ceps::getNew<Ufunc>(&solCL, CepsFunctionFlag::DofIndex));
  getFunctionDictionary()->add("cl", "OPERATOR cl_stimulation * ubar");
  CepsString sourceOptions = "cl cl UNKNOWN 0 TYPE KTERM ATTRIBUTES";
  for (CepsAttribute attr: this->getTissueAttributes())
    sourceOptions += " " + ceps::toString(attr);
  getSourceTermManager()->add(sourceOptions,false);
 
  MonodomainSolver solver(this);
  solver.solve(); 
 
  if (ceps::isProfiling())
    getProfiler()->stop("whole");  
}
 
void
CLMonodomainProblem::defineUnknowns()
{
  addUnknown(UnknownsName::vm, m_tissueAttributes);
}
 
 
void
CLMonodomainProblem::defineBoundaryConditions()
{
  MonodomainProblem::defineBoundaryConditions();
 
  CepsString stimOptions = m_parameters->getString("electrodes stimulation");
  m_functions->addFunction("cl_stimulation", stimOptions);
}