FEStaticSolver.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 "pde/solver/static/FEStaticSolver.hpp"
 
#include "pde/assemblers/AbstractAssembler.hpp"
#include "pde/problem/static/AbstractStaticPdeProblem.hpp"
 
FEStaticSolver::FEStaticSolver(AbstractStaticPdeProblem* problem) :
    AbstractStaticPdeSolver(problem),
    m_fe(ceps::runtimeCast<FiniteElements*>(problem->getSpatialDiscretization()))
{}
 
FEStaticSolver::~FEStaticSolver()
{}
 
void
FEStaticSolver::assembleAndSolve()
{
  CEPS_SAYS_NOENDL("  Assembling the linear system...");
 
  // Firstly assemble what is needed
  updateAssemblers();
 
  // Build LHS
  m_lhs->finalize();
  m_opMat->copy(*m_lhs);
  m_lhs->add(*m_bcMat, 1.0, false);
 
  // Build RHS
  m_rhs->zero();
  DVecPtr tmp = m_discretization->newDofVector();
 
  // 1. Compute source terms F and add M.F to right hand side
  if (m_hasRegSrc)
  {
    DVecPtr regSrc = m_discretization->newDofVector();
    regSrc->zero();
 
    for (auto src : m_problem->getSourceTermManager()->asVector())
      if (src->hasOption(CepsSourceTermFlag::Default) or src->hasOption(CepsSourceTermFlag::Stimulation))
        addFieldValuesTo(*src, regSrc.get(), 0.);
 
    tmp->mult(*(m_fe->getMassMatrix()), *regSrc);
    m_rhs->add(*tmp);
  }
 
  // 2. "Laplace" source terms
  if (m_hasLapSrc)
  {
    DVecPtr lapSrc = m_discretization->newDofVector();
    lapSrc->zero();
 
    for (auto src : m_problem->getSourceTermManager()->asVector())
      if (src->hasOption(CepsSourceTermFlag::Laplace))
        addFieldValuesTo(*src, lapSrc.get(), 0.);
 
    tmp->mult(*m_lapSrcMat, *lapSrc);
    m_rhs->add(*tmp);
  }
 
  // 3. The neuman and robin boundary conditions
  m_rhs->add(*m_bcVec);
 
  // Then apply Dirichlet BCs (put 1 on diagonal in the matrix, and the corresponding value into the
  // vector)
  for (auto bc : *(m_problem->getBoundaryConditionManager()->getDirichletBCs()))
    bc.second->applyAsDirichlet(m_lhs, m_rhs);
 
  m_rhs->finalize();
  m_lhs->finalize();
  
  CEPS_SAYS_INLINE(" Done");
 
  // Solve
  CEPS_SAYS_NOENDL("  Solving linear system...");
  m_linearSystem->solve(m_solution);
  CEPS_SAYS_INLINE(" Done");
}