ReferenceFE.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/finiteElements/ReferenceFE.hpp"
#include "pde/finiteElements/LagrangeBasisFunctions.hpp"
 
ReferenceFE::ReferenceFE(CepsUInt dimension, CepsFeLagrangeType type):
    m_type              (type),
    m_dim               (dimension),
    m_order             (getOrder(type)),
    m_nbDofs            (0u),
    m_nGeomVertices     (0u),
    m_ownerDeductionRule({})
{
  buildBasisFunction();
}
 
ReferenceFE::~ReferenceFE()
{}
 
CepsCellType
ReferenceFE::getCellType() const
{
  return getCategory(m_type);
}
 
CepsUInt
ReferenceFE::getPolynomialOrder() const
{
  return m_order;
}
 
CepsUInt
ReferenceFE::getNumberOfDofs() const
{
  return m_nbDofs;
}
 
CepsMathDynamic1D
ReferenceFE::computeBasisFunctions(const CepsReal3D& point)
{
  const CepsUInt &n = m_basisFunctions.size();
 
  CepsMathDynamic1D phi(n);
  for (CepsUInt i=0U; i<n; i++)
    phi(i) = computeBasisFunction(i,point);
  return phi;
}
 
CepsMathDynamic2D
ReferenceFE::computeBasisFunctionDerivatives(const CepsReal3D& point)
{
  const CepsUInt &n = m_basisFunctions.size();
 
  CepsMathDynamic2D gradPhi(n,3);
  for (CepsUInt iPhi=0U; iPhi<n; iPhi++)
  {
    const CepsReal3D &dpdx = computeBasisFunctionDerivative(iPhi, point);
    for (CepsUInt i=0; i<3; i++)
      gradPhi(iPhi,i) = dpdx[i];
  }
  return gradPhi.transpose();
}
 
const CepsVector<CepsReal3D>&
ReferenceFE::getBasisVertices() const
{
  return m_vertices;
}
 
const CepsReal3D&
ReferenceFE::getBasisVertex(CepsUInt i) const
{
  CEPS_ABORT_IF(i>=getNumberOfBasisVertices(),
    "The maximum index for basis vertices is " << getNumberOfBasisVertices()
  );
  return m_vertices[i];
}
 
CepsInt
ReferenceFE::getBasisVertexGeomCellIndex(CepsUInt i) const
{
  CEPS_ABORT_IF(i>=getNumberOfBasisVertices(),
    "The maximum index for basis vertices is " << getNumberOfBasisVertices()
  );
  return m_geomCellIndex[i];
}
 
GeomNode*
ReferenceFE::getNodeOfOwnerForBasisVertex(CepsUInt i, GeomCell* cell) const
{
  CEPS_ABORT_IF(i>=getNumberOfBasisVertices(),
    "The maximum index for basis vertices is " << getNumberOfBasisVertices()
  );
  CEPS_ABORT_IF(ceps::isNullPtr(cell),
    "Passed ptr on cell is null"
  );
 
  auto                it  = m_ownerDeductionRule[i].begin();
  GeomNode*           res = cell->getNodeAt(*it);
  CepsNodeGlobalIndex id  = res->getGlobalIndex();
 
  for (it++;it!=m_ownerDeductionRule[i].end();it++)
    if (id > cell->getNodeAt(*it)->getGlobalIndex())
    {
      res = cell->getNodeAt(*it);
      id  = res->getGlobalIndex();
    }
 
  // }
  // Never called, no quads elements implemented
  // else
  //   CEPS_ABORT("Not implemented");
 
  return res;
}
 
CepsReal3D
ReferenceFE::getTransformedBasisVertex(CepsUInt i, GeomCell* cell) const
{
  CEPS_ABORT_IF(i>=getNumberOfBasisVertices(),
    "The maximum index for basis vertices is " << getNumberOfBasisVertices()
  );
  CEPS_ABORT_IF(ceps::isNullPtr(cell),
    "Passed ptr on cell is null"
  );
  CepsMathVertex res = cell->getNodeAt(0)->getCoordinatesForEigen();
  if (m_dim>=1)
  {
    auto jac = cell->getJacobianMatrix();
    for (CepsUInt j=0;j<m_dim;j++)
      res += m_vertices[i][j]*jac.col(j);
  }
  return CepsReal3D({res[0],res[1],res[2]});
}
 
CepsUInt
ReferenceFE::getDimension() const
{
  return m_dim;
}
 
CepsUInt
ReferenceFE::getNumberOfGeomVertices() const
{
  return m_nGeomVertices;
}
 
CepsUInt
ReferenceFE::getNumberOfBasisVertices() const
{
  return m_vertices.size();
}
 
CepsUInt
ReferenceFE::getNumberOfBasisFunctions() const
{
  return m_basisFunctions.size();
}
 
CepsFeLagrangeType
ReferenceFE::getFEType() const
{
  return m_type;
}
 
void
ReferenceFE::buildBasisFunction()
{
  CepsCellType type = getCategory(m_type);
  if (type == CepsCellType::Simplex)
  {
    buildPolysBasisFunction(
      m_dim, m_order, m_vertices, m_geomCellIndex, m_basisFunctions, buildSimplexUniformPoints
    );
    m_nGeomVertices = m_dim+1;
    m_nbDofs        = m_vertices.size();
 
    if (m_dim==1) // Dimension 1 is easy, whatever the order
    {
      m_ownerDeductionRule.push_back({0});
      for (CepsUInt i=1;i<m_order;i++)
        m_ownerDeductionRule.push_back({0,1});
      m_ownerDeductionRule.push_back({1});
    }
    else 
    {
      if (m_order==1) // Order 1 is easy, whatever the dimension (nodes coincide with geom)
      {
        for (CepsUInt i=0;i<m_nbDofs;i++)
          m_ownerDeductionRule.push_back({i});
      }
      // For now, we write rules for Pk k>1 and dim>1
      else if (m_dim==2)
      {
        // 0,0 node
        m_ownerDeductionRule.push_back({0});
        // Nodes on y edge
        for (CepsUInt i=1;i<m_order;i++)
          m_ownerDeductionRule.push_back({0,2});
        // y=1 node
        m_ownerDeductionRule.push_back({2});
 
        for (CepsUInt i=1;i<m_order;i++)
        {
          // x edge
          m_ownerDeductionRule.push_back({0,1});
          // inside the triangle
          for (CepsUInt j=1;i<m_order-j;j++)
            m_ownerDeductionRule.push_back({0,2,1});
          // diagonal edge
          m_ownerDeductionRule.push_back({2,1});
        }
        // x=1 node
        m_ownerDeductionRule.push_back({1});
      }
      else if (m_dim==3)
      {
 
        // 0,0,0 node
        m_ownerDeductionRule.push_back({0});
 
        // x=0
 
        // Nodes on z edge
        for (CepsUInt i=1;i<m_order;i++)
          m_ownerDeductionRule.push_back({0,3});
        // z=1 node
        m_ownerDeductionRule.push_back({3});
 
        for (CepsUInt i=1;i<m_order;i++)
        {
          // y edge
          m_ownerDeductionRule.push_back({0,2});
          // inside the z,y triangle
          for (CepsUInt j=1;i<m_order-j;j++)
            m_ownerDeductionRule.push_back({0,3,2});
          // diagonal edge z,y
          m_ownerDeductionRule.push_back({3,2});
        }
        // y=1 node
        m_ownerDeductionRule.push_back({2});
 
        // in between slices
        for (CepsUInt i=1;i<m_order;i++)
        {
          // x edge
          m_ownerDeductionRule.push_back({0,1});
          // inside z x triangle
          for (CepsUInt j=1;j<m_order-i;j++)
            m_ownerDeductionRule.push_back({0,3,1});
          // x z edge
          m_ownerDeductionRule.push_back({3,1});
          
          // in between lines
          for (CepsUInt j=1;j<m_order-i;j++)
          {
            // y x triangle
            m_ownerDeductionRule.push_back({0,2,1});
            // in the meat of the element
            for (CepsUInt k=1;k<m_order-i-j;k++)
              m_ownerDeductionRule.push_back({0,1,2,3});
            // xyz plane
            m_ownerDeductionRule.push_back({1,2,3});
          }
 
          // xy edge
          m_ownerDeductionRule.push_back({1,2});
        }
 
        // x=1 node
        m_ownerDeductionRule.push_back({1});
      }
    }
 
  }
  // else // Quads
  // {
  //   buildPolysBasisFunction(
  //     m_dim, m_order, m_vertices, m_geomCellIndex, m_basisFunctions, buildQuadrilateralUniformPoints
  //   );
  //   m_nGeomVertices = pow(2,m_dim);
  // }
 
  return;
}
 
CepsReal
ReferenceFE::computeBasisFunction(CepsUInt iPhi, const CepsReal3D &point) 
{
  CEPS_ABORT_IF(iPhi >= getNumberOfBasisFunctions(),
    "The maximum index for basis functions is " << getNumberOfBasisFunctions()
  );
  return m_basisFunctions[iPhi].eval(point);
}
 
CepsReal3D
ReferenceFE::computeBasisFunctionDerivative(CepsUInt iPhi, const CepsReal3D &point)
{
  CEPS_ABORT_IF(iPhi >= getNumberOfBasisFunctions(),
    "The maximum index for basis functions is " << getNumberOfBasisFunctions()
  );
  return m_basisFunctions[iPhi].gradient(point);
}