FEAssembler.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/FEAssembler.hpp"
 
FEAssembler::~FEAssembler()
{
  for (auto it: m_tools)
    for (auto tool: it.second)
      ceps::destroyObject(tool);
}
 
void
FEAssembler::assemble(CepsReal t, CepsBool finalize)
{
 
  // Spatial assembly
  init();
 
  // ceps::beginSequential();
  while (not finished())
  {
    assembleForCurrentElement(t);
    next();
  }
  // ceps::endSequential();
 
  // ZeroD unknowns,
  // They are added if they were specified in the list at instanciation, or no unknowns specified,
  // if the assembler is not restricted to boundary, and if the assembler has not attribute or the
  // CepsUniversal attribute
  // Also, their contribution is not added in "integral" computation
  if (not m_boundaryOnly)
  if (this->hasAttribute(CepsUniversal) or this->getNumberOfAttributes()==0)
  for (Unknown* u: m_problem->getZeroDUnknowns())
    if (m_unknownsWL.empty() or ceps::argVector(m_unknownsWL,u)!=-1)
    {
      DegreeOfFreedom* zddof = m_fe->getDegreesOfFreedomForUnknown(u)[0];
      if (zddof->getOwner()==ceps::getRank())
        assembleForZeroDDof(zddof, t);
    }
 
  if (finalize)
  {
    if (isAssemblingMatrix())
      m_gMat->finalize();
    if (isAssemblingVector())
      m_gVec->finalize();
  }
 
  return;
}
 
FEAssembler::FEAssembler() :
  AbstractAssembler(),
  m_fe              (nullptr),
  m_boundaryOnly    (false),
  m_toParse         (nullptr),
  m_tools           ({}),
  m_bMat            ({}),
  m_bVec            ({}),
  m_dofsBlockMapping({})
{
}
 
FEAssembler::FEAssembler(FiniteElements*               fe, 
                         CepsBool                      isBoundary,
                         const CepsSet<CepsAttribute>& attributes,
                         const CepsVector<Unknown*>&   unknowns
                        ):
  AbstractAssembler(fe,attributes,unknowns),
  m_fe              (fe),
  m_boundaryOnly    (isBoundary),
  m_toParse         (nullptr),
  m_tools           ({}),
  m_bMat            ({}),
  m_bVec            ({}),
  m_dofsBlockMapping({})
{
 
  if (isBoundary)
    m_toParse = &(fe->getBoundaryFiniteElements());
  else
    m_toParse = &(fe->getFiniteElements());
  CEPS_WARNS_IF(m_toParse->empty(),
    "An assembler has no elements to parse. This can occur when the mesh has no boundary elements"
  );
 
  // Set the iterators on elements to the first to parse
  init();
 
  prepareQuadratures();
 
}
 
void
FEAssembler::prepareQuadratures()
{
 
  // No quads for now
  //for (auto cellType : {CepsCellType::Simplex,CepsCellType::Quad})
  CepsCellType cellType = CepsCellType::Simplex;
    for (CepsUInt dim = 0U; dim <= 3U; dim++)
    {
      m_tools[cellType][dim] = new IntegrationTools;
      IntegrationTools* tool = m_tools[cellType][dim];
 
      tool->ref = m_fe->getReferenceElementOfDim(cellType,dim);
 
      tool->quadrature.initializeQuadrature(tool->ref->getPolynomialOrder(),dim);
 
      const CepsUInt &n = tool->quadrature.getNbQuadPoints();
      tool->phi    .reserve(n);
      tool->gradPhi.reserve(n);
 
      for (CepsUInt i=0U; i<n; ++i)
      {
        const CepsMathVertex &xQ = tool->quadrature.getQuadPoint(i);
        const CepsReal3D     &x  = {xQ(0U),xQ(1U),xQ(2U)};
 
        tool->phi    .push_back(tool->ref->computeBasisFunctions          (x));
        tool->gradPhi.push_back(tool->ref->computeBasisFunctionDerivatives(x));
      }
 
    }
  return;
}
 
void
FEAssembler::init()
{
  if (m_toParse->empty())
    return;
 
  // The first element may not have the correct attributes
  CepsBool isUniversal = this->hasAttribute(CepsUniversal);
  m_elem = m_toParse->begin();
  if (not (isUniversal or hasOneOfAttributes((*m_elem)->getGeomObject()->getAttributes())))
    next();
 
}
 
CepsBool
FEAssembler::finished()
{
  return m_toParse->empty() or m_elem == m_toParse->end(); 
}
 
void
FEAssembler::next()
{
  CepsBool isUniversal = this->hasAttribute(CepsUniversal);
  CepsBool stop        = false;
  do
  {
    ++m_elem;
    stop = finished() or isUniversal 
    or hasOneOfAttributes((*m_elem)->getGeomObject()->getAttributes());
  }
  while (not stop); 
 
}
 
void
FEAssembler::assembleForCurrentElement(CepsReal t)
{
 
  // Fetch some information
  GeomCell* gCell = (*m_elem)->getGeomObject();
  CepsUInt  dim   = gCell->getDimension();
  auto      J     = gCell->getJacobianMatrix();
  auto      invJt = gCell->getInverseJacobianMatrix().transpose();
  auto      det   = gCell->getJacobianDeterminant();
 
  // Select the right quadrature
  IntegrationTools* itg = m_tools.at((*m_elem)->getReferenceElement()->getCellType()).at(dim);
 
  // Loop on quadrature points
  m_newBlocksNeeded = true;
 
  CepsUInt nQ = itg->quadrature.getNbQuadPoints();
  for (CepsUInt iQ = 0U; iQ < nQ; iQ++)
  {
    // Use only the necessary dimensions of the quadrature point in reference element
    auto     xQ = itg->quadrature.getQuadPoint(iQ).topLeftCorner(dim,1);
    CepsReal wQ = itg->quadrature.getWeight(iQ);
 
    // transform grad_phi using inverse jacobian
    // We use the part of grad phi which concerns the dimension of the element
    // grad_phi = (J^-1)^T * grad_phi
    auto gradPhi = invJt*itg->gradPhi[iQ].topLeftCorner(dim,(*m_elem)->getNumberOfNodes());
    auto phi     = itg->phi[iQ];
 
    // compute weight (quadrature*jacobian) and point
    CepsReal       weight = det*wQ;
    CepsMathVertex p = J*xQ + (*m_elem)->getOrigin();
    CepsReal3D     x({p[0],p[1],p[2]});
    computeBlocksOnElementAtQuadPoint(*m_elem,x,t,phi,gradPhi);
    sumBlockContribution(weight);
  }
  addBlockContribution();
 
  return;
}
 
CepsVector<Unknown*>
FEAssembler::extractUnknownsFrom(const CepsVector<Unknown*>& us, FEBase* elem) const
{
  CepsVector<Unknown*> res = {};
  for (Unknown* u: us)
  {
    if (not ceps::contains(m_unknownsWL, u) and not m_unknownsWL.empty())
      continue;
    if (not (
          u->hasUniversalAttribute() or 
          u->hasOneOfAttributes(elem->getProperties()->getAttributes())
      ))
    continue;
  
    res.push_back(u);
  }
  return res;
}
 
void
FEAssembler::setupBlocksOnElementForUnknown(FEBase* elem, Unknown* u)
{
  setupBlocksOnElementForUnknowns(elem,CepsVector<Unknown*>({u}));
  return;
}
 
void
FEAssembler::setupBlocksOnElementForUnknowns(
  FEBase* elem, 
  const CepsVector<Unknown*>& us
)
{
  // If we don't need to allocate the blocks
  // just zero them
  if (not m_newBlocksNeeded)
  {
    for (CepsUInt i = 0UL; i < m_bMat.size(); ++i)
      for (CepsUInt j = 0UL; j < m_bMat[i].size(); ++j)
        m_bMat[i][j].setZero();
    
    for (CepsUInt i = 0UL; i < m_bVec.size(); ++i)
        m_bVec[i].setZero();
 
    return;
  }
 
  CepsUInt nu    = us.size();
  CepsUInt nn    = elem->getNumberOfNodes();
  auto     attrs = elem->getGeomObject()->getAttributes();
 
  CepsVector<DegreeOfFreedom*> allDofs = m_fe->getDofsOnElementForUnknowns(elem, us);
  // allDofs was built with external loop on FE Nodes, so by taking only one per node,
  // dofs should be sorted in the same order as FE Nodes, ie the order of basis vertices.
  // \sa FiniteElements::getDofsOnElement \sa FiniteElements::buildElements
 
  m_dofsBlockMapping.clear();
  m_dofsBlockMapping.resize(nu);
  m_dofsBlock       .clear();
  m_dofsBlock       .resize(nu);
 
  for (CepsUInt i=0; i<nu; i++)
  {
    m_dofsBlockMapping[i].clear();
    m_dofsBlockMapping[i].reserve(nn);
    m_dofsBlock       [i].clear();
    m_dofsBlock       [i].reserve(nn);
    for (DegreeOfFreedom* dof:allDofs)
      if (dof->getUnknown()==us[i])
      {
        m_dofsBlock       [i].push_back(dof);
        m_dofsBlockMapping[i].push_back(dof->getGlobalIndex());
      }
  }
  
  // Prepare the submatrices
  m_bMat.clear();
  m_bMat.resize(nu);
  m_bMatSumQuads.clear();
  m_bMatSumQuads.resize(nu);
 
  m_addMatBlock.clear();
  m_addMatBlock.resize(nu);
 
  // Prepare the subvectors
  m_bVec.clear();
  m_bVec.resize(nu);
  m_bVecSumQuads.clear();
  m_bVecSumQuads.resize(nu);
 
  m_addVecBlock.clear();
  m_addVecBlock.resize(nu);
 
  CepsUInt N, M;
 
  for (CepsUInt i = 0UL; i < nu; i++)
  {
    m_bMat[i].resize(nu);
    m_bMatSumQuads[i].resize(nu);
    m_addMatBlock[i].resize(nu);
 
    for (CepsUInt j = 0UL; j < nu; j++)
    {
      N = m_dofsBlockMapping[i].size();
      M = m_dofsBlockMapping[j].size();
 
      m_bMat[i][j].setZero(N, M);
      m_bMatSumQuads[i][j].setZero(N, M);
 
      m_addMatBlock[i][j] = m_fe->getProblem()->unknownsInteract(us[i],us[j],attrs) 
                            and (N != 0) and (M != 0);
    }
 
    m_bVec[i].setZero(N);
    m_bVecSumQuads[i].setZero(N);
    m_addVecBlock[i] = true;
  }
 
  m_newBlocksNeeded = false;
}
 
void
FEAssembler::computeBlocksOnElementAtQuadPoint(
  FEBase*, CepsReal3D, CepsReal, const CepsMathDynamic1D&, 
  const CepsMathDynamic2D&)
{}
 
void 
FEAssembler::sumBlockContribution(CepsReal weight)
{
  for (CepsUInt i = 0UL; i < m_bMatSumQuads.size(); i++)
    for (CepsUInt j = 0UL; j < m_bMatSumQuads[i].size(); j++)
      m_bMatSumQuads[i][j] += weight * m_bMat[i][j];
 
  for (CepsUInt i = 0UL; i < m_bVecSumQuads.size(); i++)
    m_bVecSumQuads[i] += weight * m_bVec[i];
}
 
void
FEAssembler::addZeroDMatBlockContribution(CepsReal value, DegreeOfFreedom* a, DegreeOfFreedom* b)
{
  if (isAssemblingMatrix())
    m_gMat->addValue(m_scaleFactor * value, a->getGlobalIndex(), b->getGlobalIndex());
}
 
void
FEAssembler::addZeroDMatBlockContributions(
  const CepsVector<CepsReal>&         values,
  DegreeOfFreedom*                    a,
  const CepsVector<DegreeOfFreedom*>& b,
  CepsReal                            factor
)
{
  if (not isAssemblingMatrix())
    return;
 
  CEPS_ABORT_IF(values.size() != b.size(), "incompatible size");
  const CepsUInt& n = b.size();
  for (CepsUInt i = 0UL; i < n; ++i)
    addZeroDMatBlockContribution(factor * values[i], a, b[i]);
 
  return;
}
 
void
FEAssembler::addZeroDVecBlockContribution(CepsReal value, DegreeOfFreedom* dof)
{
  if (isAssemblingVector())
    m_gVec->addValue(m_scaleFactor * value, dof->getGlobalIndex());
}
 
void
FEAssembler::addBlockContribution()
{
  addMatrixBlockContribution();
  addVectorBlockContribution();
}
 
void
FEAssembler::addMatrixBlockContribution()
{
  if (not isAssemblingMatrix()) 
    return;
 
  for (CepsUInt i = 0UL; i < m_bMatSumQuads.size(); i++)
    for (CepsUInt j = 0UL; j < m_bMatSumQuads[i].size(); j++)
      if (m_addMatBlock[i][j]) // Add only if there is an interaction    
        m_gMat->addSubMatrix(m_scaleFactor * m_bMatSumQuads[i][j], 
          m_dofsBlockMapping[i], m_dofsBlockMapping[j]);
}
 
void
FEAssembler::addVectorBlockContribution()
{
  if (not isAssemblingVector())
    return;
 
  for (CepsUInt i = 0UL; i < m_bVecSumQuads.size(); i++)
    if (m_addVecBlock[i])
      m_gVec->addSubVector( m_scaleFactor * m_bVecSumQuads[i],
        m_dofsBlockMapping[i]);
}