FEMonodomainAssembler.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/assembler/FE/FEMonodomainAssembler.hpp"
 
FEMonodomainAssembler::FEMonodomainAssembler(
  MonodomainProblem* problem,
  FiniteElements*    fe
):
  FEHeatAssembler(fe, false, problem->getTissueAttributes(), {problem->getTMVUnknown()}),
  FECardiacParameters(problem),
  m_vm(problem->getTMVUnknown()) 
{
}
 
void
FEMonodomainAssembler::computeBlocksOnElementAtQuadPoint(
  FEBase*                  element,
  CepsReal3D               xQ,
  CepsReal                 t,
  const CepsMathDynamic1D& phi,
  const CepsMathDynamic2D& gradPhi
)
{
  setupBlocksOnElementForUnknown(element,m_vm);
 
  CepsStandardArgs args;
  args.x      = xQ;
  args.t      = t;
  args.cellId = element->getGeomObject()->getGlobalIndex();
 
  CepsMathTensor sigmaI = m_sigmaI->getValue(element->getGeomObject(),t);
  CepsMathTensor sigmaE = m_sigmaE->getValue(element->getGeomObject(),t);
  CepsMathTensor sigmaM = (sigmaI.inverse() + sigmaE.inverse()).inverse();
 
  CepsReal Am = m_Am->eval(args);
  CepsReal Cm = m_Cm->eval(args);
  
  if (getAssemblingFlag() == CepsAssemblingFlag::TimeDerivative)
  {
    m_bMat[0][0] += phi*phi.transpose();
  }
  else if (getAssemblingFlag() == CepsAssemblingFlag::SpatialOperator)
  {
    m_bMat[0][0] += 1.0/(Am*Cm)*gradPhi.transpose()*sigmaM*gradPhi;
  }
  return;
}