ExtendedBidomainSolver.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/solver/ExtendedBidomainSolver.hpp"
 
#include "cardiac/assembler/FE/FECardiacBCAssembler.hpp"
#include "cardiac/assembler/FE/FEExtendedBidomainAssembler.hpp"
#include "cardiac/problem/ExtendedBidomainProblem.hpp"
#include "geometry/selectors/AttributesSelector.hpp"
#include "pde/assemblers/FE/FENullMeanConstraintAssembler.hpp"
 
ExtendedBidomainSolver::ExtendedBidomainSolver(ExtendedBidomainProblem* problem) :
    BidomainSolver(problem)
{
  if (ceps::isValidPtr(problem->getParameters()))
    ExtendedBidomainSolver::setupWithParameters(problem->getParameters());
}
 
void
ExtendedBidomainSolver::setupWithParameters(InputParameters* params)
{
  // Scale current boundary condition
   
  auto pb = getExtendedBidomainProblem();
 
  if (pb->hasAnodalCathodalStimulation())
  {
    CEPS_SAYS("Anodal cathodal stimulation enabled");
    
    auto     bcFE        = m_fe->getBoundaryFiniteElements();
    CepsReal anodeMeas   = getDomainMeasure(bcFE, pb->getAnodeAttributes());
    CepsReal cathodeMeas = getDomainMeasure(bcFE, pb->getCathodeAttributes());
 
    BoundaryConditionManager* bcManager = pb->getBoundaryConditionManager();
 
    bcManager->getNeumannBCs()->at("ANODE")->setScaleFactor(1. / anodeMeas);
    CEPS_SAYS("Anode intensity divided by the electrode measure: " << anodeMeas);
 
    bcManager->getNeumannBCs()->at("CATHODE")->setScaleFactor(1. / cathodeMeas);
    CEPS_SAYS("Cathode intensity divided by the electrode measure: " << cathodeMeas);
  }
  
}
 
void
ExtendedBidomainSolver::initializeAssemblers()
{
  auto bepb = getExtendedBidomainProblem();
 
  if (bepb->requireNullMean())
  {
    // Lhs assembler cannot be a simple FEExtendedBidomainAssembler
    // because the ue field is defined up to a constant. Here we use the
    // actual lagragian unknown to enforce a null mean on ue
    using LhsAssembler = FENullMeanConstraintAssembler<FEExtendedBidomainAssembler>;
 
    LhsAssembler* opAsb = ceps::getNew<LhsAssembler>(bepb, m_fe);
    opAsb->addLagrangian(m_problem->getUnknown("Lagrangian for extracellular potential"));
   
    m_opAsb = opAsb;
  }
  else 
  {
    FEExtendedBidomainAssembler* opAsb = ceps::getNew<FEExtendedBidomainAssembler>(bepb,m_fe);
    
    m_opAsb = opAsb;
  }
  m_bcAsb = ceps::getNew<FECardiacBCAssembler>(bepb,CepsSet<CepsAttribute>({CepsUniversal}),m_fe);
}
 
ExtendedBidomainProblem*
ExtendedBidomainSolver::getExtendedBidomainProblem() const
{
  return ceps::runtimeCast<ExtendedBidomainProblem*>(m_problem);
}