BergotBasis.cpp

// Gmsh - Copyright (C) 1997-2018 C. Geuzaine, J.-F. Remacle
//
// See the LICENSE.txt file for license information. Please report all
// bugs and problems to the public mailing list <gmsh@onelab.info>.

#include <cmath>
#include "BergotBasis.h"
#include "MElement.h"
#include "orthogonalBasis.h"

BergotBasis::BergotBasis(int p, bool incpl) : order(p), incomplete(incpl)
{
  if(incomplete && order > 2) {
    Msg::Error("Incomplete pyramids of order %i not yet implemented", order);
  }
}

BergotBasis::~BergotBasis() {}

bool BergotBasis::validIJ(int i, int j) const
{
  if(!incomplete) return (i <= order) && (j <= order);
  if(i + j <= order) return true;
  if(i + j == order + 1) return i == 1 || j == 1;
  return false;
}

// Values of Bergot basis functions for coordinates (u,v,w) in the unit pyramid:
// f = L_i(uhat)*L_j(vhat)*(1-w)^max(i,j)*P_k^{2*max(i,j)+2,0}(what)
// with i,j = 0...order and k = 0...order-max(i,j)
// and (uhat,vhat,what) = (u/(1-w),v/(1-w),2*w-1) reduced coordinates in the
// unit cube [-1,1]^3
void BergotBasis::f(double u, double v, double w, double *val) const
{
  // Compute Legendre polynomials at uhat
  double uhat = (w == 1.) ? 0. : u / (1. - w);
  std::vector<double> uFcts(order + 1);
  LegendrePolynomials::f(order, uhat, &(uFcts[0]));

  // Compute Legendre polynomials at vhat
  double vhat = (w == 1.) ? 0. : v / (1. - w);
  std::vector<double> vFcts(order + 1);
  LegendrePolynomials::f(order, vhat, &(vFcts[0]));

  // Compute Jacobi polynomials at what
  double what = 2. * w - 1.;
  std::vector<std::vector<double> > wFcts(order + 1), wGrads(order + 1);
  for(int mIJ = 0; mIJ <= order; mIJ++) {
    int kMax = order - mIJ;
    std::vector<double> &wf = wFcts[mIJ];
    wf.resize(kMax + 1);
    JacobiPolynomials::f(kMax, 2. * mIJ + 2., 0., what, &(wf[0]));
  }

  // Recombine to find shape function values
  int index = 0;
  for(int i = 0; i <= order; i++) {
    for(int j = 0; j <= order; j++) {
      if(validIJ(i, j)) {
        int mIJ = std::max(i, j);
        double fact = pow(1. - w, mIJ);
        std::vector<double> &wf = wFcts[mIJ];
        for(int k = 0; k <= order - mIJ; k++, index++)
          val[index] = uFcts[i] * vFcts[j] * wf[k] * fact;
      }
    }
  }
}

// Derivatives of Bergot basis functions for coordinates (u,v,w) in the unit