LagrangeBasisFunctions.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/LagrangeBasisFunctions.hpp"
 
// void
// buildQuadrilateralUniformPoints(
//   CepsUInt                dim, 
//   CepsUInt                numberOfPoints, 
//   CepsVector<CepsReal3D>& points, 
//   CepsVector<CepsInt>&    geomCellIndex
// )
// {
//   CepsReal h  = 1.0 / (numberOfPoints - 1);
//   CepsUInt Nx = numberOfPoints;
//   CepsUInt Ny = (dim > 1 ? numberOfPoints : 1);
//   CepsUInt Nz = (dim > 2 ? numberOfPoints : 1);
 
//   points    .clear();
//   isGeomNode.clear();
//   points    .reserve(Nx*Ny*Nz);
//   isGeomNode.reserve(Nx*Ny*Nz);
//   for (CepsUInt i=0; i<Nx; ++i)
//     for (CepsUInt j=0; j<Ny; ++j)
//       for (CepsUInt k=0; k<Nz; ++k)
//       {
//         points.push_back(CepsReal3D({i*h, j*h, k*h}));
//         CepsInt toPush = -1;
//         if ((i==0 or i==Nx-1) and (j==0 or j==Ny-1) and (k==0 or k==Nz-1))
//           toPush = ;
//         geomCellIndex.push_back(toPush);
//       }
//   return;
// }
 
void
buildSimplexUniformPoints(
  CepsUInt                dim, 
  CepsUInt                numberOfPoints, 
  CepsVector<CepsReal3D>& points, 
  CepsVector<CepsInt>&    geomCellIndex
)
{
  CepsReal h  = 1.0 / (numberOfPoints - 1);
  CepsUInt Nx = numberOfPoints;
  CepsUInt Ny = (dim > 1 ? numberOfPoints : 1);
  CepsUInt Nz = (dim > 2 ? numberOfPoints : 1);
 
  points       .clear();
  geomCellIndex.clear();
  points       .reserve(Nx*Ny*Nz);
  geomCellIndex.reserve(Nx*Ny*Nz);
  CepsInt count = -1;
  for (CepsUInt i=0; i<Nx; ++i)
    for (CepsUInt j=0; j<Ny; ++j)
      for (CepsUInt k=0; k<Nz; ++k)
        if ((i+j+k)*h <= 1.0)  // x+y+z <= 1.0 for simplices
        {
          points.push_back(CepsReal3D({i*h, j*h, k*h}));
          CepsInt toPush = -1;
          if((i==0 or i==Nx-1) and (j==0 or j==Ny-1) and (k==0 or k==Nz-1))
          {
            toPush = geomCellIndex.size()==0 ? 0 : dim-count;
            count ++;
          }
          geomCellIndex.push_back(toPush);
        }
 
  return;
}
 
CepsVector<CepsArray3<CepsInt>>
buildBasisFunctionExponents(CepsUInt dim, CepsUInt numberOfPoints)
{
  CepsVector<CepsArray3<CepsInt>> expos = {};
 
  CepsInt Nx = numberOfPoints;  // maximum degree is numberOfPoints
  CepsInt Ny = (dim > 1 ? numberOfPoints: 1);
  CepsInt Nz = (dim > 2 ? numberOfPoints: 1);
 
  for (CepsInt i=0; i<Nx; ++i)
    for (CepsInt j=0; j<Ny; ++j)
      for (CepsInt k=0; k<Nz; ++k)
        if (i+j+k < (CepsInt)numberOfPoints)  // maximum degree is numberOfPoints
          expos.push_back(CepsArray3<CepsInt>({i,j,k}));
 
  auto compare = [] (const CepsArray3<CepsInt> &a, const CepsArray3<CepsInt> &b) -> CepsBool {
    const CepsUInt &i =
      (a[0] >= a[1] and a[0] >= a[2] ? 0U : (a[1] >= a[0] and a[1] >= a[2] ? 1U : 2U));
    const CepsUInt &j =
      (b[0] >= b[1] and b[0] >= b[2] ? 0U : (b[1] >= b[0] and b[1] >= b[2] ? 1U : 2U));
 
    const CepsInt &as = a[0] + a[1] + a[2];
    const CepsInt &bs = b[0] + b[1] + b[2];
 
    // sum are different (eg {0, 1} <= {2,0}) -> {y <= x^2}
    if (as != bs)
      return as <= bs;
    // so sum are equal here (eg {0, 2} == {1,1}) -> {x^2 ? xy}
 
    // the minimum value are different (eg {1,1} >= {0,2}) -> {xy <= x^2}
    if (a[i] != b[j])
      return a[i] <= b[j];
    // so minimum value are equal here (eg {0, 1, 1} == {1, 1, 0}) -> {xy ? yz}
 
    // only ordering the index as {x}, than {y} and than {z} (eg {2, 1, 1} <= {1, 1, 2}) -> {x^2yz
    // <= xyz^2}
    return i <= j;
  };
 
  std::sort (expos.begin(),expos.end(),compare);
 
  // only keep the number needed
  return expos;
}
 
void
buildPolysBasisFunction(
  CepsUInt                    dim,
  CepsUInt                    order,
  CepsVector<CepsReal3D>&     points,
  CepsVector<CepsInt>&        geomCellIndex,
  CepsVector<Polynomial<3U>>& polys,
  void (*pointBuilder) (CepsUInt, CepsUInt, CepsVector<CepsReal3D>&, CepsVector<CepsInt>&)
)
{
  const CepsUInt &numberOfPoints = order + 1;
 
  pointBuilder(dim,numberOfPoints,points,geomCellIndex);
  const CepsUInt &np = points.size();
 
  CepsVector<CepsArray3<CepsInt>> exponents = buildBasisFunctionExponents(dim,numberOfPoints);
  exponents.resize(np);
 
  CepsVector<Monomial<3U>> mono(np);
  for (CepsUInt i=0; i<np; ++i)
    mono[i].setExponents(exponents[i]);
 
  // mat A
  using mat_t = Eigen::Matrix<CepsReal,-1,-1>;
  mat_t A(np, np);
  for (CepsUInt i=0; i<np; ++i)
    for (CepsUInt j=0; j<np; ++j)
      A.coeffRef(i,j) = mono[j](points[i]);
 
  // second member
  mat_t I(np,np);
  I.setIdentity();
 
  // solver
  Eigen::PartialPivLU<mat_t> solver(A);
  decltype(A)                coeffs = solver.solve(I);
 
  // fill
  polys.resize(np);
  for (CepsUInt i=0; i<np; ++i)
  {
    polys[i].add(coeffs.col(i),exponents);
    polys[i].pruneZeroCoeffs();
  }
}