xref: /dports/math/xlife++/xlifepp-sources-v2.0.1-2018-05-09/src/finiteElements/triangle/NedelecTriangle.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 NedelecTriangle.cpp
  \authors D. Martin, N. Kielbasiewicz
  \since 23 nov 2004
  \date 6 aug 2012

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

#include "NedelecTriangle.hpp"
#include "../Interpolation.hpp"
#include "../integration/Quadrature.hpp"
#include "../segment/LagrangeSegment.hpp"
#include "utils.h"
#include <iostream>

namespace xlifepp
{

//! NedelecTriangle constructor for Nedelec reference elements
NedelecTriangle::NedelecTriangle(const Interpolation* interp_p)
  : RefTriangle(interp_p)
{
  name_ += "_Nedelec";
  mapType = _covariantPiolaMap;
  dofCompatibility=_signDofCompatibility;
  dimShapeFunction=2;
}

NedelecTriangle::~NedelecTriangle() {}

/* return the dof permutation when triangle vertices are permuted
   i,j,k are the vertex numbers in local numbering (1,2 or 3)
   when (i,j,k) = (1,2,3)  no dof permutation (return a void vector!)
   Example with a 2 order element

     v2                     v3
      | \                    | \
    3 *  * 2               4 *  * 5
      |     \                |     \            perm =[ 2 1 6 5 4 3 7 8]
    4 *  *7   * 1          3 *  *7   * 6
      |    *8   \            |    *8   \
     v3---*---*--v1         v2---*---*--v1
          5   6                  2   1
    original numbering     (i,j,k) = (2,1,3)

    No permutation of internal dofs !!!  to be checked
*/
std::vector<number_t> NedelecTriangle::dofsMap(const number_t& i, const number_t& j, const number_t& k) const
{
  if(i==1 && j==2) return std::vector<number_t>();  //no permutation
  number_t d = interpolation_p->numtype;
  number_t ni=3*d+1;
  std::vector<number_t> perm(d*(d+2));
  std::vector<number_t>::iterator p=perm.begin();
  if(i==1 && j==3)  // (i,j,d) = (1,3,2) : d=1 -> [3 2 1]  d=2 -> [6 5 4 3 2 1 7 8]
    {
      number_t n = 3*d;
      for(number_t j=0; j<3*d; j++) *p++=n--;
      n=3*d+1;
      for(number_t j=0; j<d*(d-1); j++) *p++=ni++;
      return perm;
    }
  if(i==2)
    {
      if(j==1) // (i,j,d) = (2,1,3) : d=1 -> [1 3 2]  d=2 -> [2 1 6 5 4 3 7 8]
        {
          number_t n = d;
          for(number_t j=0; j<d; j++) *p++=n--;
          n=3*d;
          for(number_t j=0; j<d; j++) *p++=n--;
          n=2*d;
          for(number_t j=0; j<d; j++) *p++=n--;
          for(number_t j=0; j<d*(d-1); j++) *p++=ni++;
        }
      else // (i,j,d) = (2,3,1) : d=1 -> [2 3 1]  d=2 -> [3 4 5 6 1 2 7 8]
        {
          number_t n = d+1;
          for(number_t j=0; j<d; j++) *p++=n++;
          n=2*d+1;
          for(number_t j=0; j<d; j++) *p++=n++;
          n=1;
          for(number_t j=0; j<d; j++) *p++=n++;
          for(number_t j=0; j<d*(d-1); j++) *p++=ni++;

        }
      return perm;
    }
  //i=3
  if(j==1) // (i,j,d) = (3,1,2) : d=1 -> [3 1 2]  d=2 -> [5 6 1 2 3 4 7 8]
    {
      number_t n = 2*d+1;
      for(number_t j=0; j<d; j++) *p++=n++;
      n=1;
      for(number_t j=0; j<d; j++) *p++=n++;
      n=d+1;
      for(number_t j=0; j<d; j++) *p++=n++;
      for(number_t j=0; j<d*(d-1); j++) *p++=ni++;
    }
  else // (i,j,d) = (3,2,1) : d=1 -> [2 1 3]  d=2 -> [4 3 2 1 6 5 7 8]
    {
      number_t n = 2*d;
      for(number_t j=0; j<d; j++) *p++=n--;
      n=d;
      for(number_t j=0; j<d; j++) *p++=n--;
      n=3*d;
      for(number_t j=0; j<d; j++) *p++=n--;
      for(number_t j=0; j<d*(d-1); j++) *p++=ni++;
    }
  return perm;
}


//! Nedelec standard P1 triangle Reference Element
NedelecFirstTriangleP1::NedelecFirstTriangleP1(const Interpolation* interp_p)
  : NedelecTriangle(interp_p)
{
  name_ += "_first family_1";
  interpolationData();  // build element interpolation data
  sideNumbering();      // local numbering of d.o.f's on edges
  pointCoordinates();   // "virtual coordinates" of moment dofs
  //sideRefElement();     // Reference Element on edges
  maxDegree = 1;
}

NedelecFirstTriangleP1::~NedelecFirstTriangleP1() {}

//! interpolationData defines Reference Element interpolation data
void NedelecFirstTriangleP1::interpolationData()
{
  trace_p->push("NedelecTriangle::interpolationData");
  nbDofsInSides_ = nbDofs_ = 3;
  // Creating Degrees Of Freedom of reference element, in Nedelec standard elements all D.o.F on edges are sharable
  refDofs.reserve(nbDofs_);
  for(number_t s = 1; s <= 3; s++)
    refDofs.push_back(new RefDof(true, _onEdge, s, 1, 2, 1, 0, 2, 0, _crossnProjection, "int_e u.tq"));
  shapeValues.resize(*this, 2); // allocate storage for shape functions
  nbPts_ = nbDofs_;
  trace_p->pop();
}

/* builds virtual coordinates of moment dofs of order k
     for side dofs : {1/2}
*/
void NedelecFirstTriangleP1::pointCoordinates()
{
  std::vector<RefDof*>::iterator it_rd=refDofs.begin();
  (*it_rd++)->coords(0.5,0.5);
  (*it_rd++)->coords(0.,0.5);
  (*it_rd++)->coords(0.5,0.);
}

//! sideNumbering defines Reference Element local numbering on sides
void NedelecFirstTriangleP1::sideNumbering()
{
  trace_p->push("NedelecTriangle::sideNumbering");
  /*
     Local numbering of Nedelec P1 first family

         *
         | \
         2   1
         |     \
         *---3---*
         .. k=1
   */
  sideDofNumbers_.resize(3);
  for(number_t s = 0; s < 3; s++) sideDofNumbers_[s].push_back(s + 1);
  trace_p->pop();
}

void NedelecFirstTriangleP1::computeShapeValues(std::vector<real_t>::const_iterator it_pt, ShapeValues& shv, const bool withDeriv) const
{
  /*
     Shape functions    : w_1 = (-y, x) , w_2 = (-y, -1+x) , w_3 = (1-y, x)
     Degrees of Freedom :
   */
  real_t x1 = *it_pt, x2 = *(it_pt + 1);
  std::vector<real_t>::iterator it_w(shv.w.begin());
  *it_w++ = -x2;
  *it_w++ = x1; // w_1
  *it_w++ = -x2;
  *it_w++ = x1 - 1.; // w_2
  *it_w++ = 1. - x2;
  *it_w   = x1; // w_3

  if(withDeriv)
    {
      std::vector<std::vector<real_t> >::iterator it_dwdx1(shv.dw.begin()), it_dwdx2(it_dwdx1 + 1);
      std::vector<real_t>::iterator it_dwdx1i((*it_dwdx1).begin()), it_dwdx2i((*it_dwdx2).begin());
      *it_dwdx1i++ = 0.; *it_dwdx2i++ = -1.;
      *it_dwdx1i++ = 1.; *it_dwdx2i++ = 0.; // grad w_1
      *it_dwdx1i++ = 0.; *it_dwdx2i++ = -1.;
      *it_dwdx1i++ = 1.; *it_dwdx2i++ = 0.; // grad w_2
      *it_dwdx1i++ = 0.; *it_dwdx2i++ = -1.;
      *it_dwdx1i   = 1.; *it_dwdx2i   = 0.; // grad w_3
    }
}

//=====================================================================================
/*! Nedelec first family order k on triangle

        *             *
        | \           | \
        2   1       3 *   * 2
        |     \       |     \           dof's side numbering
        *---3---*   4 *  *7   * 1
          k=1         |    *8   \
                      *---*---*---*
                          5   6
                           k=2
        2
        | \
        2   1        edge orientation in dof's integral computation on edge
        |     \
        3---3---1         e1 : 1->2,  e2 : 2->3   e3 : 3->1
           ->
*/
//=====================================================================================
NedelecFirstTrianglePk::NedelecFirstTrianglePk(const Interpolation* interp_p)
  : NedelecTriangle(interp_p)
{
  name_ += "_first family_"+tostring(interp_p->numtype);
  interpolationData();  // build element interpolation data
  sideNumbering();      // local numbering of points on edges
  pointCoordinates();   // "virtual coordinates" of moment dofs
  //sideRefElement();     // Reference Element on edges
}

NedelecFirstTrianglePk::~NedelecFirstTrianglePk() {}

//! interp defines Reference Element interpolation data
void NedelecFirstTrianglePk::interpolationData()
{
  trace_p->push("NedelecFirstTrianglePk::interpolationData");

  number_t k = interpolation_p->numtype;
  nbDofs_= k*(k+2);
  nbDofsInSides_=3*k;
  nbInternalDofs_=k*(k-1);

  /* Creating Degrees Of Freedom of reference element
     in Nedelec first family order k elements, D.o.Fs are defined as moment D.o.Fs
        edge dofs     : v-> int_e v.t q,  q in P_(k-1)[e]     k dofs by edge e
        internal dofs : v-> int_t v.q,    q in P_(k-2)[t]^2   k(k-1) dofs  only for k>1
     numbering of D.oFs are given by the basis functions order of dual spaces defining D.o.Fs
     NOTE : moment D.o.Fs can be reinterpreted as punctual D.o.Fs using a well adapted quadrature formulae
            we do not use this representation here
  */
  refDofs.reserve(nbDofs_);
  number_t nbds=nbDofsInSides_/3;
  for(number_t e = 1; e <= 3; ++e)              //edge dofs (rank in dual polynomial basis is used to identify dof)
    for(number_t i=0; i<nbds; ++i)
      refDofs.push_back(new RefDof(true, _onEdge, e, i+1, 2, 1, 0, 2, 0, _crossnProjection, "int_e u.t q"));
  for(number_t i=0; i<nbInternalDofs_; ++i) //element dofs (order = i is used to identify dof)
    refDofs.push_back(new RefDof(false, _onElement, 0, i+1, 2, 2, 0, 2, 0, _noProjection, "int_k u.q"));

  // compute shape functions as formal polynomials
  computeShapeFunctions();
  buildPolynomialTree();   //construct tree representation
  maxDegree = Wk.degree();

  // allocate storage for shape functions
  shapeValues.resize(*this, 2);

  trace_p->pop();
}

//! compute shape functions as polynomials using general algorithm
void NedelecFirstTrianglePk::computeShapeFunctions()
{
  number_t k = interpolation_p->numtype;
  PolynomialsBasis::iterator itp;

  //compute Matrix Lij=l(pi,qj) for all pi in Vk, for all qj in P_(k-1)[e], edge D.o.Fs :  v-> int_e v.t q, q in P_(k-1)[e]
  PolynomialsBasis Vk(_Rk,2,k);       //Nedelec first family polynomials space (Rk)
  LagrangeStdSegment segPk(findInterpolation(_Lagrange, _standard, k-1, H1));
  segPk.computeShapeFunctions();
  PolynomialsBasis& Q=segPk.Wk;       //create P_(k-1)[x]
  PolynomialsBasis::iterator itq;
  dimen_t d=Vk.degree()+Q.degree();
  QuadRule bestQR = Quadrature::bestQuadRule(_segment,d);
  Quadrature* quad=findQuadrature(_segment, bestQR, d);  //quadrature rule on segment of order d (integrates exactly x^d)
  number_t nbq= quad->numberOfPoints();
  Matrix<real_t> L(nbDofs_,nbDofs_);
  number_t j=1;
  real_t sq2=std::sqrt(2)/2.;
  std::vector<real_t>::const_iterator itpt, itw; //for quadrature rule
  for(itp=Vk.begin(); itp!=Vk.end(); ++itp, ++j)
    {
      number_t i=1;
      for(number_t e = 1; e <= 3; ++e)
        {
          Vector<real_t> t(2,-1.);
          if(e==1) {t[1]=1; t*=sq2;}        // on edge 1 : t = (-1,1)/sqrt(2)
          if(e==2) t[0]=0;                  // on edge 2 : t = (0,-1)
          if(e==3) {t[0]=1; t[1]=0;}        // on edge 3 : t = (1,0)
          Polynomial pn=dot(*itp,t);
          for(itq=Q.begin(); itq!=Q.end(); ++itq, ++i)
            {
              L(i,j)=0;
              itpt=quad->point();
              itw=quad->weight();
              Polynomial& pq=(*itq)[0];
              for(number_t q=0; q<nbq; ++q, ++itpt, ++itw)
                {
                  real_t s=*itpt;
                  if(e==1)  L(i,j)+= pn(1-s,s) * pq(s)** itw * 2* sq2;
                  if(e==2)  L(i,j)+= pn(0,1-s) * pq(s)** itw ;
                  if(e==3)  L(i,j)+= pn(s,0)   * pq(s)** itw ;
                }
            }
        }
    }
  if(k>1) //compute Matrix Lij=l(pi,qj) for all pi in Vk, qj in P_(k-2)[t]^2, element  D.o.Fs :  v-> int_t v.q,  q in P_(k-2)[t]^2
    {
      PolynomialsBasis R(PolynomialBasis(_Pk,2,k-2),2); //P_(k-2)[t]^2
      PolynomialsBasis::iterator itr;
      dimen_t d=Vk.degree()+R.degree();
      bestQR = Quadrature::bestQuadRule(_triangle,d);
      quad=findQuadrature(_triangle,bestQR, d);  //quadrature rule on triangle of order d
      nbq= quad->numberOfPoints();
      j=1;
      for(itp=Vk.begin(); itp!=Vk.end(); ++itp,++j)
        {
          number_t i=nbDofsInSides_+1;

          for(itr=R.begin(); itr!=R.end(); ++itr,++i)
            {
              L(i,j)=0.;
              itpt=quad->point();
              itw=quad->weight();
              for(number_t q=0; q<nbq; ++q, itpt+=2, ++itw)
                {
                  real_t x=*itpt, y=*(itpt+1);
                  L(i,j)+= dot(eval(*itp,x,y),eval(*itr,x,y))** itw ;
                }
            }
        }
    }
  //compute coefficients of shape functions by solving linear systems
  Matrix<real_t> sf(nbDofs_,_idMatrix);
  Matrix<real_t> Lo=L;
  real_t eps=theZeroThreshold;
  number_t r;
  bool ok=gaussMultipleSolver(L, sf, nbDofs_, eps, r);
  if(!ok)
    {
      where("NedelecFirstTrianglePk::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=2;
  Wk.dimVec=2;
  Wk.name="RT"+tostring(k);
  for(number_t i=0; i<nbDofs_; i++)
    {
      for(itv=v.begin(); itv!=v.end(); ++itv, ++itm) *itv=*itm;
      Wk.push_back(combine(Vk,v));
    }
  //Wk.clean(0.000001);
  dWk.resize(2);
  dWk[0]=dx(Wk);
  dWk[1]=dy(Wk);
}

/* builds virtual coordinates of moment dofs of order k
     for side dofs : {1/k+1, 2/k+1, ..., k/k+1}
     for internal dofs : coordinates of Pk-2 on triangle (1/k+1, 1/k+1), (k-1/k+1, 1/k+1), (1/k+1, k-1/k+1),
                         2 internal dofs (i,i+1) located at the same place
*/
void NedelecFirstTrianglePk::pointCoordinates()
{
  std::vector<RefDof*>::iterator it_rd=refDofs.begin();
  number_t nbds=nbDofsInSides_/3, k= interpolation_p->numtype, kp=k+1;
  for(number_t e = 1; e <= 3; ++e)              //edge dofs
    for(number_t i=1; i<=nbds; ++i, ++it_rd)
      {
        real_t a = real_t(i)/kp, b=1.-a;
        switch(e)
          {
            case 1: (*it_rd)->coords(b, a); break;
            case 2: (*it_rd)->coords(0, b); break;
            default : (*it_rd)->coords(a, 0);
          }
      }
  for(number_t j=1; j<=k-1; ++j)   //internal dofs
    {
      real_t y = real_t(j)/kp;
      for(number_t i=1; i<=k-j; ++i)
        {
          real_t x = real_t(i)/kp;
          (*it_rd++)->coords(x, y);
          (*it_rd++)->coords(x, y);
        }
    }
}

//dof's side numbering
void NedelecFirstTrianglePk::sideNumbering()
{
  trace_p->push("NedelecFirstTrianglePk::sideNumbering");
  number_t k = interpolation_p->numtype;
  number_t n=1;
  sideDofNumbers_.resize(3,std::vector<number_t>(k,0));
  for(number_t s = 0; s < 3; ++s)          // k dofs by side (edge)
    for(number_t i=0; i<k; ++i, ++n)
      sideDofNumbers_[s][i] = n;
  trace_p->pop();
}

//tool for global dofs construction
//return new index from a given one n from local edge vertex numbers i,j (k not used)
number_t NedelecFirstTrianglePk::sideDofsMap(const number_t& n, const number_t& i, const number_t& j, const number_t& k) const
{
  number_t d=interpolation_p->numtype;
  if(d==1 || i<j) return n;   //no reverse
  if(n<=2) return 3-n;        //"extremal" side dofs
  return d+3-n;               //"internal" side dofs
}

//compute shape functions using polynomial representation
void NedelecFirstTrianglePk::computeShapeValues(std::vector<real_t>::const_iterator it_pt, ShapeValues& shv, const bool withDeriv) const
{
  computeShapeValuesFromShapeFunctions(it_pt, shv, withDeriv);
}

//=====================================================================================
/*! Nedelec second family order k on triangle

                      |\
        *             4  3
        2 4           |    \
        |   \         5      2
        5     1       |   10   \           dof's side numbering
        *-3--6-*      6  11 12   1
          k=1         |            \
                      *--7---8---9---
                           k=2
        2
        | \
        2   1        edge orientation in dof's integral computation on edge
        |     \
        3---3---1         e1 : 1->2,  e2 : 1->3   e3 : 3->1
           ->
*/
//=====================================================================================
NedelecSecondTrianglePk::NedelecSecondTrianglePk(const Interpolation* interp_p)
  : NedelecTriangle(interp_p)
{
  name_ += "_second family_"+tostring(interp_p->numtype);
  interpolationData();  // build element interpolation data
  sideNumbering();      // local numbering of points on edges
  pointCoordinates();   // set virtual dof coordinates
}

NedelecSecondTrianglePk::~NedelecSecondTrianglePk() {}

//! interp defines Reference Element interpolation data
void NedelecSecondTrianglePk::interpolationData()
{
  trace_p->push("NedelecSecondTrianglePk::interpolationData");

  number_t k = interpolation_p->numtype;
  nbDofs_= (k+1)*(k+2);
  nbDofsInSides_=3*(k+1);
  nbInternalDofs_=k*k-1;

  /* Creating Degrees Of Freedom of reference element
     in Nedelec second family order k elements, D.o.Fs are defined as moment D.o.Fs
        edge dofs     : v-> int_e v.t q,  q in P_k[e]     (k+1) dofs by edge e
        internal dofs : v-> int_t v.q,    q in D_(k-1)[t]  k*k-1 dofs  only for k>1
     numbering of D.oFs are given by the basis functions order of dual spaces defining D.o.Fs
     NOTE : moment D.o.Fs can be reinterpreted as punctual D.o.Fs using a well adapted quadrature formulae
            we do not use this representation here
  */
  refDofs.reserve(nbDofs_);
  number_t nbds=nbDofsInSides_/3;
  for(number_t e = 1; e <= 3; ++e)              //edge dofs (rank in dual polynoms basis is used to identify dof)
    for(number_t i=0; i<nbds; ++i)
      refDofs.push_back(new RefDof(true, _onEdge, e, i+1, 2, 1, 0, 2, 0, _crossnProjection, "int_e u.t q"));
  for(number_t i=0; i<nbInternalDofs_; ++i) //element dofs (order = i is used to identify dof)
    refDofs.push_back(new RefDof(false, _onElement, 0, i+1, 2, 2, 0, 2, 0, _noProjection, "int_k u.q"));

  // compute shape functions as formal polynomials
  computeShapeFunctions();
  buildPolynomialTree();   //construct tree representation

  // allocate storage for shape functions
  shapeValues.resize(*this, 2);

  trace_p->pop();
}

//! compute shape functions as polynomials using general algorithm
void NedelecSecondTrianglePk::computeShapeFunctions()
{
  number_t k = interpolation_p->numtype;
  PolynomialsBasis::iterator itp;

  //compute Matrix Lij=l(pi,qj) for all pi in Vk, for all qj in P_(k-1)[e], edge D.o.Fs :  v-> int_e v.t q, q in P_(k-1)[e]
  PolynomialsBasis Vk(PolynomialBasis(_Pk,2,k),2);       //(Pk)^2
  LagrangeStdSegment segPk(findInterpolation(_Lagrange, _standard, k, H1));
  segPk.computeShapeFunctions();
  PolynomialsBasis& Q=segPk.Wk;       //create P_(k)[x]
  PolynomialsBasis::iterator itq;
  dimen_t d=Vk.degree()+Q.degree();
  QuadRule bestQR = Quadrature::bestQuadRule(_segment,d);
  Quadrature* quad=findQuadrature(_segment, bestQR, d);  //quadrature rule on segment of order d (integrates exactly x^d)
  number_t nbq= quad->numberOfPoints();
  Matrix<real_t> L(nbDofs_,nbDofs_);
  number_t j=1;
  real_t sq2=std::sqrt(2)/2.;
  std::vector<real_t>::const_iterator itpt, itw; //for quadrature rule
  for(itp=Vk.begin(); itp!=Vk.end(); ++itp, ++j)
    {
      number_t i=1;
      for(number_t e = 1; e <= 3; ++e)
        {
          Vector<real_t> t(2,-1.);
          if(e==1) {t[1]=1; t*=sq2;}        // on edge 1 : t = (-1,1)/sqrt(2)
          if(e==2) t[0]=0;                  // on edge 2 : t = (0,-1)
          if(e==3) {t[0]=1; t[1]=0;}        // on edge 3 : t = (1,0)
          Polynomial pn=dot(*itp,t);
          for(itq=Q.begin(); itq!=Q.end(); ++itq, ++i)
            {
              L(i,j)=0;
              itpt=quad->point();
              itw=quad->weight();
              Polynomial& pq=(*itq)[0];
              for(number_t q=0; q<nbq; ++q, ++itpt, ++itw)
                {
                  real_t s=*itpt;
                  if(e==1)  L(i,j)+= pn(1-s,s) * pq(s)** itw * 2* sq2;
                  if(e==2)  L(i,j)+= pn(0,1-s) * pq(s)** itw ;
                  if(e==3)  L(i,j)+= pn(s,0)   * pq(s)** itw ;
                }
            }
        }
    }
  if(k>1) //compute Matrix Lij=l(pi,qj) for all pi in Vk, qj in D_(k-1)[t], element  D.o.Fs :  v-> int_t v.q,  q in D_(k-1)[t]
    {
      PolynomialsBasis R(_Dk, 2, k-1); //D_(k-1)[t]
      PolynomialsBasis::iterator itr;
      dimen_t d=Vk.degree()+R.degree();
      bestQR = Quadrature::bestQuadRule(_triangle,d);
      quad=findQuadrature(_triangle, bestQR, d);  //quadrature rule on triangle of order d
      nbq= quad->numberOfPoints();
      j=1;
      for(itp=Vk.begin(); itp!=Vk.end(); ++itp,++j)
        {
          number_t i=nbDofsInSides_+1;

          for(itr=R.begin(); itr!=R.end(); ++itr,++i)
            {
              L(i,j)=0.;
              itpt=quad->point();
              itw=quad->weight();
              for(number_t q=0; q<nbq; ++q, itpt+=2, ++itw)
                {
                  real_t x=*itpt, y=*(itpt+1);
                  L(i,j)+= dot(eval(*itp,x,y),eval(*itr,x,y))** itw ;
                }
            }
        }
    }
  //compute coefficients of shape functions by solving linear systems
  Matrix<real_t> sf(nbDofs_,_idMatrix);
  Matrix<real_t> Lo=L;
  real_t eps=theZeroThreshold;
  number_t r;
  bool ok=gaussMultipleSolver(L, sf, nbDofs_, eps, r);
  if(!ok)
    {
      where("NedelecSecondTrianglePk::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=2;
  Wk.dimVec=2;
  Wk.name="RT"+tostring(k);
  for(number_t i=0; i<nbDofs_; i++)
    {
      for(itv=v.begin(); itv!=v.end(); ++itv, ++itm) *itv=*itm;
      Wk.push_back(combine(Vk,v));
    }
  Wk.clean(0.000001);
  dWk.resize(2);
  dWk[0]=dx(Wk);
  dWk[1]=dy(Wk);
}

/* builds virtual coordinates of moment dofs of order k
     for side dofs : {1/k+1, 2/k+1, ..., k/k+1}
     for internal dofs : set to (1/3,1/3), to be improved in future
*/
void NedelecSecondTrianglePk::pointCoordinates()
{
  std::vector<RefDof*>::iterator it_rd=refDofs.begin();
  number_t nbds=nbDofsInSides_/3, k= interpolation_p->numtype, kp=k+2;
  for(number_t e = 1; e <= 3; ++e)              //edge dofs
    for(number_t i=1; i<=nbds; ++i, ++it_rd)
      {
        real_t a = real_t(i)/kp, b=1.-a;
        switch(e)
          {
            case 1: (*it_rd)->coords(b, a); break;
            case 2: (*it_rd)->coords(0, b); break;
            default : (*it_rd)->coords(a, 0);
          }
      }

  for(number_t i=0; i<nbInternalDofs_; ++i) //For the moment all are set to (1/3,1/3)
    (*it_rd++)->coords(1./3, 1./3);
}

//dof's side numbering
void NedelecSecondTrianglePk::sideNumbering()
{
  trace_p->push("NedelecSecondTrianglePk::sideNumbering");
  number_t k = interpolation_p->numtype;
  number_t n=1;
  sideDofNumbers_.resize(3,std::vector<number_t>(k,0));
  for(number_t s = 0; s < 3; ++s)          // k dofs by side (edge)
    for(number_t i=0; i<k; ++i, ++n)
      sideDofNumbers_[s][i] = n;
  trace_p->pop();
}

//tool for global dofs construction
//return new index from a given one n from local side numbers i,j (k not used)
number_t NedelecSecondTrianglePk::sideDofsMap(const number_t& n, const number_t& i, const number_t& j, const number_t& k) const
{
  if(i<j) return n;                      //no reverse
  if(n<=2) return 3-n;                   //"extremal" side dofs
  return (interpolation_p->numtype+3)-n; //"internal" side dofs
}

//compute shape functions using polynomial representation
void NedelecSecondTrianglePk::computeShapeValues(std::vector<real_t>::const_iterator it_pt, ShapeValues& shv, const bool withDeriv) const
{
  computeShapeValuesFromShapeFunctions(it_pt, shv, withDeriv);
}

} // end of namespace xlifepp