finiteElements/segment/LagrangeSegment.cpp

/*
XLiFE++ is an extended library of finite elements written in C++
    Copyright (C) 2014  Lunéville, Eric; Kielbasiewicz, Nicolas; Lafranche, Yvon; Nguyen, Manh-Ha; Chambeyron, Colin

    This program 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.
    This program 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 this program.  If not, see <http://www.gnu.org/licenses/>.
 */

/*!
  \file LagrangeSegment.cpp
  \authors D. Martin, N. Kielbasiewicz
  \since 17 dec 2002
  \date 6 aug 2012

  \brief Implementation of xlifepp::LagrangeSegment class members and related functions
 */

#include "LagrangeSegment.hpp"
#include "../Interpolation.hpp"
#include "utils.h"
#include "mathsResources.h"

#include <iostream>

namespace xlifepp
{

//! LagrangeSegment constructor for Lagrange reference elements
LagrangeSegment::LagrangeSegment(const Interpolation* interp_p)
   : RefSegment(interp_p)
{
  name_ += "_Lagrange";
  trace_p->push("LagrangeSegment::LagrangeSegment (" + name_ + ")");
  // build element interpolation data
  interpolationData();
  // local numbering of points on sides
  sideNumbering();
  // Reference Element on edges
  sideRefElement();
  maxDegree = interpolation_p->numtype;

  trace_p->pop();
}

LagrangeSegment::~LagrangeSegment() {}

//! interpolationData defines Lagrange Reference Element interpolation data
void LagrangeSegment::interpolationData()
{
  trace_p->push("LagrangeSegment::interpolationData");
  number_t kd = interpolation_p->numtype;

  switch (kd)
  {
    case 0: nbDofs_ = nbInternalDofs_ = 1; break;
    case 1: nbDofs_ = nbDofsOnVertices_ = 2; break;
    default: // Lagrange segment of degree > 1
      nbDofs_ = kd + 1;
      nbInternalDofs_ = kd - 1;
      nbDofsOnVertices_ = 2;
      break;
  }

  // creating D.o.F of reference element
  // in Lagrange standard elements only vertices are sharable
  refDofs.reserve(nbDofs_);
  lagrangeRefDofs(1, nbDofsOnVertices_, nbInternalDofs_);
  //lagrangeRefDofs(1, nbDofsOnVertices_, nbDofs_ - nbInternalDofs_, 1);  //old method

  // allocate storage for shape functions
  shapeValues.resize(*this, 1);
  nbPts_ = nbDofs_;
  trace_p->pop();
}

//! create standard Lagrange numbering on sides
void LagrangeSegment::sideNumbering()
{
  number_t nbSides = geomRefElem_p->nbSides();
  sideDofNumbers_.resize(nbSides);
  sideDofNumbers_[0].resize(1);
  sideDofNumbers_[1].resize(1);

  sideDofNumbers_[0][0] = geomRefElem_p->sideVertexNumber(1, 1);
  sideDofNumbers_[1][0] = geomRefElem_p->sideVertexNumber(1, 2);
}

//! return nodes numbers of first order segment elements when splitting current element
std::vector<std::vector<number_t> > LagrangeSegment::splitP1() const {
  std::vector<std::vector<number_t> > splitnum;
  std::vector<number_t> pair(2,0);
  number_t k=interpolation_p->numtype;
  switch (k) {
    case 0:
    case 1:
      pair[0]=1; pair[1]=2;
      splitnum.push_back(pair);
      return splitnum;
    case 2:
      pair[0]=1; pair[1]=3;
      splitnum.push_back(pair);
      pair[0]=3; pair[1]=2;
      splitnum.push_back(pair);
      return splitnum;
    case 3:
      pair[0]=1; pair[1]=3;
      splitnum.push_back(pair);
      pair[0]=3; pair[1]=4;
      splitnum.push_back(pair);
      pair[0]=4; pair[1]=2;
      splitnum.push_back(pair);
      return splitnum;
    case 4:
      pair[0]=1; pair[1]=3;
      splitnum.push_back(pair);
      pair[0]=3; pair[1]=4;
      splitnum.push_back(pair);
      pair[0]=4; pair[1]=5;
      splitnum.push_back(pair);
      pair[0]=5; pair[1]=2;
      splitnum.push_back(pair);
      return splitnum;
    case 5:
      pair[0]=1; pair[1]=3;
      splitnum.push_back(pair);
      pair[0]=3; pair[1]=4;
      splitnum.push_back(pair);
      pair[0]=4; pair[1]=5;
      splitnum.push_back(pair);
      pair[0]=5; pair[1]=6;
      splitnum.push_back(pair);
      pair[0]=6; pair[1]=2;
      splitnum.push_back(pair);
      return splitnum;
    default:
      pair[0]=1; pair[1]=3;
      splitnum.push_back(pair);
      for (number_t i=0; i<k-2; i++) {
        pair[0]=i+3; pair[1]=i+4;
        splitnum.push_back(pair);
      }
      pair[0]=k+1; pair[1]=2;
      splitnum.push_back(pair);
      return splitnum;
  }
  return splitnum;
}

//! return nodes numbers of first order segment elements when splitting current element
splitvec_t LagrangeSegment::splitO1() const
{
  splitvec_t splitnum;
  std::vector<std::vector<number_t> > buf;
  buf=this->splitP1();
  for (number_t i=0; i< buf.size(); i++)
  {
    splitnum.push_back(std::make_pair(_segment, buf[i]));
  }
  return splitnum;
}

//! output of Lagrange D.o.f's numbers on underlying P1 D.o.f's
void LagrangeSegment::outputAsP1(std::ofstream& os, const int refNum, const std::vector<number_t>& ranks) const
{
  number_t kd = interpolation_p->numtype;
  switch (kd)
  {
    case 0: noSuchFunction("outputAsP1"); break;
    case 1: simplexVertexOutput(os, refNum, ranks[0], ranks[1]); break;
    default:
      simplexVertexOutput(os, refNum, ranks[0], ranks[2]);
      for (number_t i = 2; i < kd; i++)
      {
        simplexVertexOutput(os, refNum, ranks[i], ranks[i + 1]);
      }
      simplexVertexOutput(os, refNum, ranks[kd], ranks[1]);
  }
}

//! Lagrange segment Reference Element for Pk interpolation at equidistant nodes
LagrangeStdSegment::LagrangeStdSegment(const Interpolation* interp_p)
   : LagrangeSegment(interp_p)
{
  name_ += "_" + tostring(interp_p->numtype);
  // local coordinates of points supporting D.o.F
  pointCoordinates();
  // build O1 splitting scheme
  splitO1Scheme = splitO1();
}

/*!
  pointCoordinates defines Lagrange Reference Element point coordinates

  Local numbering of Lagrange Pk 1D elements
  \verbatim
  --1-- 2---1  2---3---1  2---4---3---1  2---5---4---3---1  2---6---5---4---3---1
   k=0   k=1      k=2          k=3              k=4                  k=5
  \endverbatim
 */
void LagrangeStdSegment::pointCoordinates()
{
  number_t kd = interpolation_p->numtype;
  std::vector<RefDof*>::iterator it_rd(refDofs.begin());
  real_t x=1.;
  real_t unsurk=1./kd;

  switch (kd)
  {
    case 0:  // Lagrange [Std] segment of order kd=0
      (*it_rd)->coords(.5);
      break;
    case 1:  // Lagrange [Std] segment of order kd=1
      (*it_rd++)->coords(1.);
      (*it_rd++)->coords(0.);
      break;
    case 2:  // Lagrange [Std] segment of order kd=2
      (*it_rd++)->coords(1.);
      (*it_rd++)->coords(0.);
      (*it_rd++)->coords(0.5);
      break;
    default: // Lagrange [Std] segment of order kd > 2
      (*it_rd++)->coords(1.);
      (*it_rd++)->coords(0.);
      while (it_rd != refDofs.end())
      {
        x -= unsurk;
        (*it_rd++)->coords(x);
      }
      break;
  }
}

/*!
  shapeFunctions defines Lagrange Reference Element shape functions

  Shape functions (Lagrange polynomials at \f$n_0=1, n_1=0, n_2=1-1/n, n_3=1-2/n, \ldots n_n=1/n\f$)
   with \f$w_i(n_j)=\delta_{ij} \; (0 <= i,j <= n)\f$, are given by

  \f$\begin{array}{rl}
                     w_0 & \displaystyle = \frac{n^{n-1}}{(n-1)!} x \prod_{j=1,n-1} (x-j/n) \\
                         & \displaystyle = \frac{x}{(n-1)!} \prod_{j=1,n-1} (nx-j) \\
                     w_1 & \displaystyle = (-1)^n n^{n-1}C_n^n (x-1)\frac{x-1/n}{1}\ldots\frac{x-(n-1)/n}{n-1} \\
                         & \displaystyle = (-1)^n \frac{n^{n-1}}{(n-1)!} (x-1) \prod_{j=1,n-1} (x-j/n) \\
                         & \displaystyle = (-1)^n \frac{x-1}{(n-1)!} \prod_{j=1,n-1} (nx-j) \\
                     w_2 & \displaystyle = (-1) n^{n-1} C_n^1 x(x-1)\frac{(x-1/n)}{1}\ldots\frac{(x-(j-1)/n)}{(j-1)}\frac{(x-j)/n)}{j}*\frac{(x-(j+1)/n)}{(j+1)}\ldots\frac{(x-(n-2)/n)}{(n-2)} \\
                         & \displaystyle = (-1) \frac{n^{n-1}}{(n-1)!} \frac{n!}{((n-2)2!)} \prod_{j=0,n \& j!= n-1} (x-j/n) \\
                         & \displaystyle = (-1) \frac{nx(x-1)}{((n-2)2!)} \prod_{j=1,n-2} (nx-j) \\
    \forall i=3,n \; w_i & \displaystyle = (-1)^{i-1} n^{n-1} C_n^{i-1} x(x-1)\frac{(x-1/n)}{1}\ldots\frac{(x-(k-1)/n)}{(k-1)}\frac{(x-(k+1)/n)}{(k+1)}\ldots\frac{(x-(n-1)/n)}{(n-1)} \\
                         & \displaystyle = (-1)^{i-1} \frac{n^{n-1}}{(n-1)!} \frac{n!}{((n-i)i!)} \prod_{j=1,n-1 \& j!=k} (x-j/n) \\
                         & \displaystyle = (-1)^{i-1} \frac{n*x(x-1)}{((n-i)i!)} \prod_{j=1,n-1 \& j!=k} (n*x-j)
  \end{array}\f$

  where \f$k=n+1-i\f$ and \f$\displaystyle C_n^j=\frac{n!}{(n-j)! j!}\f$ are Pascal's triangle binomial coefficients.

  Examples :
  - n=1
    - \f$ w_0 = 1*1^0 x = x\f$
    - \f$ w_1 = -1*1^(x-1) = 1-x\f$
  - n=2
    - \f$ w_0 = 1*2^1/1 \;x (x-1/2) = x *(2x-1)\f$
    - \f$ w_1 = 1*2^1/1 \;(x-1)(x-1/2) = (x-1)*(2x-1)\f$
    - \f$ w_2 = -2*2^1/1 \;x(x-1) = -4x(x-1)\f$
  - n=3
    - \f$ w_0 = 1*3^2/2 \; x(x-1/3)(x-2/3) = 1/2 \; x(3x-1)(3x-2)\f$
    - \f$ w_1 = -1*3^2/2 \; (x-1)(x-1/3)(x-2/3) = -1/2 \; (x-1)(3x-1)(3x-2)\f$
    - \f$ w_2 = -3*3^2/2 \; x(x-1)(x-1/3) = -9/2 \; x(x-1)(3x-1)\f$
    - \f$ w_3 = 3*3^2/2 \; x(x-1)(x-2/3)= 9/2 \; x(x-1)(3x-2)\f$
 */
void LagrangeStdSegment::computeShapeValues(std::vector<real_t>::const_iterator it_pt, ShapeValues& shv, const bool withDeriv) const
{
  // Lagrange standard shape function computation for any order kd >= 0
  // kd : interpolation order
  int kd = interpolation_p->numtype, kp1_i;
  real_t x = *it_pt, kdx = kd * x, tmp;
  std::vector<real_t>::iterator it_w(shv.w.begin());
  std::vector<std::vector<real_t> >::iterator it_dw(shv.dw.begin());
  std::vector<real_t>::iterator it_dwi((*it_dw).begin());
  switch (kd)
  {
    case 0:  // w_0=1
      *it_w = 1.;
      if (withDeriv) {*it_dwi = 0.;}
      break;
    case 1:  // w_0=x , w_1=1-x
      *it_w++ = x;
      *it_w = 1. - x;
      if (withDeriv)
      {
        *it_dwi++ = 1.;
        *it_dwi = -1.;
      }
      break;
    case 2: // w_0=x*(2x-1) , w_1=(x-1)*(2x-1) , w_2=-4x(x-1)
      *it_w++ = x * (kdx - 1.);
      *it_w++ = (x - 1.) * (kdx - 1.);
      *it_w = -4.*x * (x - 1.);
      if (withDeriv)
      {
        *it_dwi++ = 4 * x - 1.;
        *it_dwi++ = 4 * x - 3.;
        *it_dwi = -8 * x + 4.;
      }
      break;
    default:
      //  shape function general computation
      //
      // NOTE : Below we have  x_0=1, x_1=0, xl=1-(l-1)/k=(k+1-l)/k for l=2..k
      //
      // w_0=(x-0)/(1-0) * Prod_{l=2, k} (x-xl)/(1-xl) with w_0(1)=1.
      //    =x * Prod_{l=2..k} (kx -(k+1-l))/(k -(k+1-l))
      *it_w = x;
      for (int j = 1; j < kd; j++) {*it_w *= (kdx - j) / (kd - j);}
      // w_1=(x-1)/(0-1) * Prod_{l=2..k} (x-xl)/(-xl) with w_1(0)=1.
      it_w++;
      *it_w = 1. - x;
      for (int j = 1; j < kd; j++) {*it_w *= -(kdx - j) / j;}
      // shape functions w_i, i=2..k
      kp1_i = kd - 1;
      for (it_w = shv.w.begin() + 2; it_w != shv.w.end(); kp1_i--, it_w++)
      {
        // w_i=(x-1/xi-1) * (x-0/xi-0)  * Prod_{l=2..k, l!=i} (x-xl)/(xi-xl)
        *it_w = 1.;
        for (int l = 0; l < kp1_i; l++) {*it_w *= (kdx - l) / (kp1_i - l);}
        for (int l = kp1_i + 1; l <= kd; l++) {*it_w *= (kdx - l) / (kp1_i - l);}
      }

      if (withDeriv)
      {
        //    First derivative of shape functions
        //
        // derivative of shape functions w_i, i=0, k+1
        // compute product of values (x-xl / xi-xl) for l=0..k+1, l!=i,j
        // then multiply by (1 / xi-xj) then make sum on j's

        // derivative of shape function dw_0
        *it_dwi = 0.;
        for (int j = 0; j < kd; j++)
        {
          tmp = 1.;
          for (int l = 0; l < j; l++) {tmp *= (kdx - l) / (kd - l);}
          for (int l = j + 1; l < kd; l++) {tmp *= (kdx - l) / (kd - l);}
          *it_dwi += tmp * kd / (kd - j);
        }

        // derivative of shape function dw_1
        it_dwi++;
        *it_dwi = 0.;
        for (int j = 1; j <= kd; j++)
        {
          tmp = 1.;
          for (int l = 1; l < j; l++) {tmp *= -(kdx - l) / l;}
          for (int l = j + 1; l <= kd; l++) {tmp *= -(kdx - l) / l;}
          *it_dwi -= tmp * kd / j;
        }

        // derivative of shape functions dw_i (0 <i <= kd)
        kp1_i = kd - 1;
        for (it_dwi = (*it_dw).begin() + 2; it_dwi != (*it_dw).end(); kp1_i--, it_dwi++)
        {
          *it_dwi = 0.;
          for (int j = 0; j < kp1_i; j++)
          {
            tmp = 1.;
            for (int l = 0; l < j; l++) {tmp *= (kdx - l) / (kp1_i - l);}
            for (int l = j + 1; l < kp1_i; l++) {tmp *= (kdx - l) / (kp1_i - l);}
            for (int l = kp1_i + 1; l <= kd; l++) {tmp *= (kdx - l) / (kp1_i - l);}
            *it_dwi += tmp * kd / (kp1_i - j);
          }
          for (int j = kp1_i + 1; j <= kd; j++)
          {
            tmp = 1.;
            for (int l = 0; l < kp1_i; l++) {tmp *= (kdx - l) / (kp1_i - l);}
            for (int l = kp1_i + 1; l < j; l++) {tmp *= (kdx - l) / (kp1_i - l);}
            for (int l = j + 1; l <= kd; l++) {tmp *= (kdx - l) / (kp1_i - l);}
            *it_dwi += tmp * kd / (kp1_i - j);
          }
        }
      }
  }
}

//! compute shape functions as polynomials using general algorithm
// not used by default
void LagrangeStdSegment::computeShapeFunctions()
{
  //compute Matrix Lij=pj(ni) for all pi in Vk, ni nodes
  number_t k=interpolation_p->numtype;
  Matrix<real_t> L(nbDofs_,nbDofs_);
  PolynomialBasis Pk(_Pk,1,k);             //Pk[X] space equipped with canonical basis
  std::vector<RefDof*>::iterator itd=refDofs.begin();
  for (number_t i = 1; i <= nbDofs_; ++i, ++itd)   //loop on nodes
  {
    std::vector<real_t>::const_iterator itpt=(*itd)->coords();
    PolynomialBasis::iterator itp = Pk.begin();
    for (number_t j=1; j<=nbDofs_; ++j, ++itp) L(i,j)=itp->eval(*itpt);  //loop on monomials
  }
  //compute coefficients of shape functions by solving linear systems
  Matrix<real_t> sf(nbDofs_,_idMatrix);
  real_t eps=theZeroThreshold;
  number_t r;
  bool ok=gaussMultipleSolver(L, sf, nbDofs_, eps, r);
  if (!ok)
  {
    where("LagrangeTrianglePk::computeShapeFunctions()");
    error("mat_noinvert");
  }
  sf.roundToZero();
  Matrix<real_t>::iterator itm=sf.begin();
  Vector<real_t> v(nbDofs_);
  Vector<real_t>::iterator itv;
  Wk.dimVar=1;   Wk.dimVec=1;
  Wk.name="P"+tostring(k);
  for (number_t i=0; i<nbDofs_; i++)
  {
    for (itv=v.begin(); itv!=v.end(); ++itv, ++itm) *itv=*itm;
    Wk.add(combine(Pk,v));
  }
  //Wk.clean(0.000001);
  dWk.resize(2);
  dWk[0]=dx(Wk); dWk[1]=dy(Wk);
}

//! Lagrange segment Reference Element for Pk interpolation at Gauss-Lobatto nodes
LagrangeGLSegment::LagrangeGLSegment(const Interpolation* interp_p)
   : LagrangeSegment(interp_p)
{
  name_ += "_Gauss-Lobatto_" + tostring(interp_p->numtype);
  // local coordinates of points supporting D.o.F
  pointCoordinates();
  // build O1 splitting scheme
  splitO1Scheme = splitO1();
}

//! pointCoordinates defines Lagrange Reference Element point coordinates
void LagrangeGLSegment::pointCoordinates()
{
  number_t kd(interpolation_p->numtype);
  std::vector<RefDof*>::iterator it_rd(refDofs.begin());
  std::vector<real_t> coord1d;

  switch (kd)
  {
    case 0:  // Lagrange Gauss-Lobatto segment of order kd=0 (unavailable)
      break;
    case 1:  // Lagrange [Gauss-Lobatto=Std] segment of order kd=1
      (*it_rd++)->coords(1.);
      (*it_rd++)->coords(0.);
      break;
    case 2:  // Lagrange [Gauss-Lobatto=Std] segment of order kd=2
      (*it_rd++)->coords(1.);
      (*it_rd++)->coords(0.);
      (*it_rd++)->coords(0.5);
      break;
    default: // Lagrange Gauss-Lobatto segment of order kd > 2
      coord1d.resize((nbDofs_ + 1) / 2);
      gaussLobattoPoints(nbDofs_, coord1d);
      (*it_rd++)->coords(1.);
      (*it_rd++)->coords(0.);
      // GaussLobattoPoints(n) returns the first (n+1)/2 non-negative points
      // of Gauss-Lobatto computed on [-1,1] in ascending order,
      // including middle point of abscissa 0 for n odd and end point of abscissa 1.
      // build point 3 up to (n+1)/2, for n > 3
      for (std::vector<real_t>::reverse_iterator c1d = coord1d.rbegin() + 1; c1d != coord1d.rend(); c1d++) {(*it_rd++)->coords((1. + *c1d) / 2);}
      // build points (n+3)/2 up to n
      for (std::vector<real_t>::const_iterator c1d = coord1d.begin() + nbDofs_ % 2; c1d != coord1d.end() - 1; c1d++) {(*it_rd++)->coords((1. - *c1d) / 2);}
      break;
  }
}

//! shapeFunctions defines Lagrange Reference Element shape functions
void LagrangeGLSegment::computeShapeValues(std::vector<real_t>::const_iterator it_pt, ShapeValues& shv, const bool withDeriv) const
{
  // Lagrange Gauss-Lobatto shape function computation for order kd > 0
  real_t x(*it_pt), xi, xj, xl, tmp;
  std::vector<RefDof*>::const_iterator it_rdb(refDofs.begin()), it_rde(refDofs.end()), it_rdi, it_rdl, it_rdj;
  std::vector<real_t>::iterator it_w(shv.w.begin());
  //  shape function general computation
  //
  // NOTE : Below we have  x_0=1, x_1=0 and xl (l=2..k) are the other Gauss-Lobatto
  // points in ]0,1[ ordered by decreasing abscissa.
  //
  // w_0=(x-0)/(1-0) * Prod_{l=2..k} (x-xl)/(1-xl) with w_0(1)=1.
  *it_w = x;
  for (it_rdl = it_rdb + 2; it_rdl != it_rde; it_rdl++) {xl = *((*it_rdl)->coords()); *it_w *= (x - xl) / (1. - xl);}
  // w_1=(x-1)/(0-1) * Prod_{l=2..k} (x-xl)/(-xl) with w_1(0)=1.
  it_w++;
  *it_w = 1. - x;
  for (it_rdl = it_rdb + 2; it_rdl != it_rde; it_rdl++) {xl = *((*it_rdl)->coords()); *it_w *= -(x - xl) / xl;}
  // shape functions w_i, i=2..k
  it_rdi = it_rdb + 2;
  for (it_w = shv.w.begin() + 2; it_w != shv.w.end(); it_rdi++, it_w++)
  {
    // w_i=(x-1)(xi-1) * (x-0/xi-0) * Prod_{l=2..k, l!=i} (x-xl)/(xi-xl)
    xi = *((*it_rdi)->coords());
    // compute (x-1)(xi-1) * (x-0/xi-0)
    *it_w = ((x - 1.) / (xi - 1.)) * (x / xi);
    // then multiply by Prod_{l=2..k, l!=i} (x-xl)/(xi-xl)
    for (it_rdl = it_rdb + 2; it_rdl != it_rdi; it_rdl++) {xl = *((*it_rdl)->coords()); *it_w *= (x - xl) / (xi - xl);}
    for (it_rdl = it_rdi + 1; it_rdl != it_rde; it_rdl++) {xl = *((*it_rdl)->coords()); *it_w *= (x - xl) / (xi - xl);}
  }

  if (withDeriv)
  {
    //    Derivative of shape functions
    //
    // derivative of shape functions w_i, i=0..k
    // compute product of values (x-xl / xi-xl) for l=0..k, l!=i,j
    // multiply by (1 / xi-xj) then sum on j's
    std::vector<std::vector<real_t> >::iterator it_dw(shv.dw.begin());
    std::vector<real_t>::iterator it_dwi;
    // loop on i=0..kd (all shape functions)
    it_rdi = it_rdb;

    for (std::vector<real_t>::iterator it_dwi = (*it_dw).begin(); it_dwi != (*it_dw).end(); it_rdi++, it_dwi++)
    {
      *it_dwi = 0.;
      xi = *((*it_rdi)->coords());
      // loop on j=0..i-1
      for (it_rdj = it_rdb; it_rdj != it_rdi; it_rdj++)
      {
        xj = *((*it_rdj)->coords());
        tmp = 1.;
        // loops on l=0..kd , l!=j , l!=i (j <i)
        for (it_rdl = it_rdb; it_rdl != it_rdj; it_rdl++) {xl = *((*it_rdl)->coords()); tmp *= (x - xl) / (xi - xl);}
        for (it_rdl = it_rdj + 1; it_rdl != it_rdi; it_rdl++) {xl = *((*it_rdl)->coords()); tmp *= (x - xl) / (xi - xl);}
        for (it_rdl = it_rdi + 1; it_rdl != it_rde; it_rdl++) {xl = *((*it_rdl)->coords()); tmp *= (x - xl) / (xi - xl);}
        *it_dwi += tmp / (xi - xj);
      }

      // loop on j=i+1..kd
      for (it_rdj = it_rdi + 1; it_rdj != it_rde; it_rdj++)
      {
        xj = *((*it_rdj)->coords());
        tmp = 1.;
        // loops on l=0..kd , l!=i , l!=j (j > i)
        for (it_rdl = it_rdb; it_rdl != it_rdi; it_rdl++) {xl = *((*it_rdl)->coords()); tmp *= (x - xl) / (xi - xl);}
        for (it_rdl = it_rdi + 1; it_rdl != it_rdj; it_rdl++) {xl = *((*it_rdl)->coords()); tmp *= (x - xl) / (xi - xl);}
        for (it_rdl = it_rdj + 1; it_rdl != it_rde; it_rdl++) {xl = *((*it_rdl)->coords()); tmp *= (x - xl) / (xi - xl);}
        *it_dwi += tmp / (xi - xj);
      }
    }
  }
}

} // end of namespace xlifepp