FEDivKGradAssembler.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/assemblers/FE/FEDivKGradAssembler.hpp"
 
FEDivKGradAssembler::FEDivKGradAssembler():
  FEAssembler(),
  m_cstScalars({}),
  m_fctScalars({}),
  m_cstTensors({}),
  m_fctTensors({}),
  m_isTimeDependant(false)
{
}
 
FEDivKGradAssembler::FEDivKGradAssembler(
  FiniteElements*               fe, 
  CepsBool                      isBoundary,
  const CepsSet<CepsAttribute>& attributes,
  const CepsVector<Unknown*>&   unknowns
):
  FEAssembler(fe,isBoundary,attributes,unknowns),
  m_cstScalars({}),
  m_fctScalars({}),
  m_cstTensors({}),
  m_fctTensors({}),
  m_isTimeDependant(false)
{
}
 
FEDivKGradAssembler::~FEDivKGradAssembler()
{
}
 
void
FEDivKGradAssembler::setKForUnknown(Unknown* u, CepsMathScalar k)
{
  m_cstScalars[u->getIdentifier()] = k;
}
 
void
FEDivKGradAssembler::setKForUnknown(Unknown* u, SAFunc<CepsMathScalar>* k)
{
  m_fctScalars[u->getIdentifier()] = k;
 
  if (k->hasOption(CepsFunctionFlag::Time))
    m_isTimeDependant = true;
}
 
void
FEDivKGradAssembler::setKForUnknown(Unknown* u, CepsMathTensor k)
{
  m_cstTensors[u->getIdentifier()] = k;
}
 
void
FEDivKGradAssembler::setKForUnknown(Unknown* u, SAFunc<CepsMathTensor>* k)
{
  m_fctTensors[u->getIdentifier()] = k;
  
  if (k->hasOption(CepsFunctionFlag::Time))
    m_isTimeDependant = true;
}
 
CepsBool
FEDivKGradAssembler::isChangingBetweenTimes(CepsReal t1, CepsReal t2) const
{
  return m_isTimeDependant;
 
}
 
CepsVector<CepsMathTensor>
FEDivKGradAssembler::getDiffusionTensorsForUnknowns(
  const CepsVector<Unknown*>& us, 
  FEBase*                     elem,
  CepsReal3D                  x, 
  CepsReal                    t
) const
{
  CepsVector<CepsMathTensor> res;
  res.reserve(us.size());
  for (Unknown* u: us)
  {
    if (u->isOnLocation(CepsLocationFlag::ZeroD))
    {
      res.push_back({});
      continue;
    }
    
    // Loop through all 4 kinds of diffusion coefficients
 
    CepsUnknownIndex uId = u->getIdentifier();
    // Constant tensors, not much to do
    if (m_cstTensors.contains(uId))
      res.push_back(m_cstTensors.at(uId));
    // Constant scalars, multipliy by identity matrix
    else if (m_cstScalars.contains(uId))
    {
      CepsMathDynamic2D id = Eigen::MatrixXd::Identity(3,3);
      res.push_back(m_cstScalars.at(uId)*id);
    }
    else // Functions, we need to build the arguments
    {
      CepsStandardArgs args;
      args.x      = x;
      args.t      = t;
      args.cellId = elem->getGeomObject()->getGlobalIndex();
      args.attr   = elem->getGeomObject()->getAttributes();
 
      if (m_fctTensors.contains(uId))
        res.push_back(m_fctTensors.at(u->getIdentifier())->operator()(args));
      // Scalar function: multiply by identity matrix
      else if (m_fctScalars.contains(uId))
      {
        CepsMathDynamic2D id = Eigen::MatrixXd::Identity(3,3);
        res.push_back(m_fctScalars.at(u->getIdentifier())->operator()(args)*id);
      }
      else 
      {
        CEPS_ABORT(
          "Tensor for unknown " << std::quoted(u->getName()) << " was not set." << std::endl <<
          "   Use setKForUnknown(...) before assembly."
        );
      }
    }
 
  }
 
  return res;
}
 
void
FEDivKGradAssembler::computeBlocksOnElementAtQuadPoint(
  FEBase*                  elem,
  CepsReal3D               xQ,
  CepsReal                 t,
  const CepsMathDynamic1D& phi,
  const CepsMathDynamic2D& gradPhi
)
{
  // Add Gradphi*GradphiT for each spatial unknown
  auto us       = extractUnknownsFrom(m_problem->getUnknowns(), elem);
  auto tensors  = getDiffusionTensorsForUnknowns(us,elem,xQ,t);
 
  setupBlocksOnElementForUnknowns(elem,us);
 
  for (CepsUInt i=0; i<us.size(); i++)
    if (us[i]->isSpatial())
      m_bMat[i][i] += gradPhi.transpose() * tensors[i] * gradPhi;
 
  return;
}