add MetricBasis

Numeric/CMakeLists.txt

@@ -16,6 +16,7 @@ set(SRC
 	legendrePolynomials.cpp
     bezierBasis.cpp
     JacobianBasis.cpp
+    MetricBasis.cpp
     pointsGenerators.cpp
   ElementType.cpp
   GaussIntegration.cpp

Numeric/JacobianBasis.cpp

@@ -605,6 +605,17 @@ void JacobianBasis::getMetricMinAndGradients(const fullMatrix<double> &nodesXYZ,
 
 }
 
+void JacobianBasis::interpolate(const fullVector<double> &jacobian,
+                                const fullMatrix<double> &uvw,
+                                fullMatrix<double> &result) const
+{
+  fullMatrix<double> bezM(jacobian.size(), 1);
+  fullVector<double> bez;
+  bez.setAsProxy(bezM, 0);
+  lag2Bez(jacobian, bez);
+  bezier->interpolate(bezM, uvw, result);
+}
+
 fullMatrix<double> JacobianBasis::generateJacMonomialsPyramid(int order)
 {
   int nbMonomials = (order+3)*((order+3)+1)*(2*(order+3)+1)/6 - 5;

Numeric/JacobianBasis.h

@@ -72,22 +72,22 @@ class JacobianBasis {
   inline void getSignedJacobian(const fullMatrix<double> &nodesXYZ, fullVector<double> &jacobian) const {
     getSignedJacobianGeneral(numJacNodes,gradShapeMatX,gradShapeMatY,gradShapeMatZ,nodesXYZ,jacobian);
   }
-  inline void getSignedJacobianFast(const fullMatrix<double> &nodesXYZ, fullVector<double> &jacobian) const {
-    getSignedJacobianGeneral(numJacNodesFast,gradShapeMatXFast,gradShapeMatYFast,gradShapeMatZFast,nodesXYZ,jacobian);
-  }
   inline void getSignedJacobian(const fullMatrix<double> &nodesX, const fullMatrix<double> &nodesY,
                                 const fullMatrix<double> &nodesZ, fullMatrix<double> &jacobian) const {
     getSignedJacobianGeneral(numJacNodes,gradShapeMatX,gradShapeMatY,
                              gradShapeMatZ,nodesX,nodesY,nodesZ,jacobian);
   }
+  inline void getScaledJacobian(const fullMatrix<double> &nodesXYZ, fullVector<double> &jacobian) const {
+    getScaledJacobianGeneral(numJacNodes,gradShapeMatX,gradShapeMatY,gradShapeMatZ,nodesXYZ,jacobian);
+  }
+  inline void getSignedJacobianFast(const fullMatrix<double> &nodesXYZ, fullVector<double> &jacobian) const {
+    getSignedJacobianGeneral(numJacNodesFast,gradShapeMatXFast,gradShapeMatYFast,gradShapeMatZFast,nodesXYZ,jacobian);
+  }
   inline void getSignedJacobianFast(const fullMatrix<double> &nodesX, const fullMatrix<double> &nodesY,
                                     const fullMatrix<double> &nodesZ, fullMatrix<double> &jacobian) const {
     getSignedJacobianGeneral(numJacNodesFast,gradShapeMatXFast,gradShapeMatYFast,
                              gradShapeMatZFast,nodesX,nodesY,nodesZ,jacobian);
   }
-  inline void getScaledJacobian(const fullMatrix<double> &nodesXYZ, fullVector<double> &jacobian) const {
-    getScaledJacobianGeneral(numJacNodes,gradShapeMatX,gradShapeMatY,gradShapeMatZ,nodesXYZ,jacobian);
-  }
   inline void getScaledJacobianFast(const fullMatrix<double> &nodesXYZ, fullVector<double> &jacobian) const {
     getScaledJacobianGeneral(numJacNodesFast,gradShapeMatXFast,gradShapeMatYFast,gradShapeMatZFast,nodesXYZ,jacobian);
   }
@@ -104,6 +104,10 @@ class JacobianBasis {
   inline void subdivideBezierCoeff(const fullVector<double> &bez, fullVector<double> &result) const {
     bezier->subDivisor.mult(bez,result);
   }
+  //
+  void interpolate(const fullVector<double> &jacobian,
+                   const fullMatrix<double> &uvw,
+                   fullMatrix<double> &result) const;
 
   // Jacobian basis order and pyramidal basis
   void getGradientsFromNodes(const fullMatrix<double> &nodes,

Numeric/MetricBasis.cpp

+// Gmsh - Copyright (C) 1997-2013 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@geuz.org>.
+
+#include "MetricBasis.h"
+//#include "GmshDefines.h"
+#include "BasisFactory.h"
+#include "pointsGenerators.h"
+#include <cmath>
+
+MetricCoefficient::MetricCoefficient(MElement *el) : _element(el)
+{
+  _jacobian = BasisFactory::getJacobianBasis(el->getTypeForMSH());
+
+  int nJac = _jacobian->getNumJacNodes();
+  int nMapping = _jacobian->getNumMapNodes();
+  fullMatrix<double> nodesXYZ(nMapping, 3);
+  el->getNodesCoord(nodesXYZ);
+
+  switch (el->getDim()) {
+  case 0 :
+    return;
+
+  case 1 :
+  {
+    Msg::Fatal("not implemented");
+  }
+  break;
+
+  case 2 :
+  {
+    Msg::Fatal("not implemented");
+  }
+  break;
+
+  case 3 :
+  {
+    fullMatrix<double> dxyzdX(nJac,3), dxyzdY(nJac,3), dxyzdZ(nJac,3);
+    _jacobian->getGradientsFromNodes(nodesXYZ, &dxyzdX, &dxyzdY, &dxyzdZ);
+
+    _coefficientsLag.resize(nJac, 7);
+    for (int i = 0; i < nJac; i++) {
+      const double &dxdX = dxyzdX(i,0), &dydX = dxyzdX(i,1), &dzdX = dxyzdX(i,2);
+      const double &dxdY = dxyzdY(i,0), &dydY = dxyzdY(i,1), &dzdY = dxyzdY(i,2);
+      const double &dxdZ = dxyzdZ(i,0), &dydZ = dxyzdZ(i,1), &dzdZ = dxyzdZ(i,2);
+      const double dvxdX = dxdX*dxdX + dydX*dydX + dzdX*dzdX;
+      const double dvxdY = dxdY*dxdY + dydY*dydY + dzdY*dzdY;
+      const double dvxdZ = dxdZ*dxdZ + dydZ*dydZ + dzdZ*dzdZ;
+      _coefficientsLag(i, 0) = (dvxdX + dvxdY + dvxdZ) / 3;
+      _coefficientsLag(i, 1) = dvxdX - _coefficientsLag(i, 0);
+      _coefficientsLag(i, 2) = dvxdY - _coefficientsLag(i, 0);
+      _coefficientsLag(i, 3) = dvxdZ - _coefficientsLag(i, 0);
+      _coefficientsLag(i, 4) = dxdX*dxdY + dydX*dydY + dzdX*dzdY;
+      _coefficientsLag(i, 5) = dxdZ*dxdY + dydZ*dydY + dzdZ*dzdY;
+      _coefficientsLag(i, 6) = dxdX*dxdZ + dydX*dydZ + dzdX*dzdZ;
+    }
+  }
+  break;
+  }
+}
+
+void MetricCoefficient::getCoefficients(fullMatrix<double> &coeff, bool bezier)
+{
+  fullMatrix<double> *ref;
+  switch (_coefficientsLag.size2()) {
+  case 1:
+    if (!bezier) coeff = _coefficientsLag;
+    else {
+      if (_coefficientsBez.size2() != 1) _computeBezCoeff();
+      coeff = _coefficientsBez;
+    }
+    break;
+
+  case 3:
+    if (!bezier) ref = &_coefficientsLag;
+    else {
+      if (_coefficientsBez.size2() != 3) _computeBezCoeff();
+      ref = &_coefficientsBez;
+    }
+    coeff.resize(ref->size1(), 2);
+    for (int i = 0; i < ref->size1(); ++i) {
+      double tmp = pow((*ref)(i, 1), 2);
+      tmp += pow((*ref)(i, 2), 2);
+      tmp = std::sqrt(tmp);
+      coeff(i, 0) = (*ref)(i, 0) - tmp;
+      coeff(i, 1) = (*ref)(i, 0) + tmp;
+    }
+    break;
+
+  case 7:
+    if (!bezier) ref = &_coefficientsLag;
+    else {
+      if (_coefficientsBez.size2() != 7) _computeBezCoeff();
+      ref = &_coefficientsBez;
+    }
+    coeff.resize(ref->size1(), 2);
+    for (int i = 0; i < ref->size1(); ++i) {
+      double tmp = pow((*ref)(i, 1), 2);
+      tmp += pow((*ref)(i, 2), 2);
+      tmp += pow((*ref)(i, 3), 2);
+      tmp += 2 * pow((*ref)(i, 4), 2);
+      tmp += 2 * pow((*ref)(i, 5), 2);
+      tmp += 2 * pow((*ref)(i, 6), 2);
+      tmp = std::sqrt(tmp);
+      double factor = std::sqrt(6)/3;
+      coeff(i, 0) = (*ref)(i, 0) - factor * tmp;
+      coeff(i, 1) = (*ref)(i, 0) + factor * tmp;
+    }
+    break;
+
+  default:
+    Msg::Error("Wrong number of functions for metric: %d",
+               _coefficientsLag.size2());
+  }
+}
+
+static int nChoosek(int n, int k)
+{
+  if (n < k || k < 0) {
+    Msg::Error("Wrong argument for combination.");
+    return 1;
+  }
+
+  if (k > n/2) k = n-k;
+  if (k == 1)
+    return n;
+  if (k == 0)
+    return 1;
+
+  int c = 1;
+  for (int i = 1; i <= k; i++, n--) (c *= n) /= i;
+  return c;
+}
+
+void MetricCoefficient::interpolate(const double *uvw, double *minmaxQ)
+{
+  if (minmaxQ == NULL) {
+    Msg::Error("Cannot write solution of interpolation");
+    return;
+  }
+
+  int order = _jacobian->bezier->getOrder();
+
+  int dimSimplex;
+  fullMatrix<double> exponents;
+  double bezuvw[3];
+  switch (_element->getType()) {
+  case TYPE_PYR:
+    bezuvw[0] = .5 * (1 + uvw[0]);
+    bezuvw[1] = .5 * (1 + uvw[1]);
+    bezuvw[2] = uvw[2];
+    _interpolateBezierPyramid(uvw, minmaxQ);
+    return;
+
+  case TYPE_HEX:
+    bezuvw[0] = .5 * (1 + uvw[0]);
+    bezuvw[1] = .5 * (1 + uvw[1]);
+    bezuvw[2] = .5 * (1 + uvw[2]);
+    dimSimplex = 0;
+    exponents = gmshGenerateMonomialsHexahedron(order);
+    break;
+
+  case TYPE_TET:
+    bezuvw[0] = uvw[0];
+    bezuvw[1] = uvw[1];
+    bezuvw[2] = uvw[2];
+    dimSimplex = 3;
+    exponents = gmshGenerateMonomialsTetrahedron(order);
+    break;
+
+  case TYPE_PRI:
+    bezuvw[0] = uvw[0];
+    bezuvw[1] = uvw[1];
+    bezuvw[2] = .5 * (1 + uvw[2]);
+    dimSimplex = 2;
+    exponents = gmshGenerateMonomialsPrism(order);
+    break;
+  }
+
+  int numCoeff = exponents.size1();
+  int dim = exponents.size2();
+
+  if (_coefficientsBez.size2() == 0) _computeBezCoeff();
+
+  double *terms = new double[_coefficientsBez.size2()];
+  for (int t = 0; t < _coefficientsBez.size2(); ++t) {
+    terms[t] = 0;
+    for (int i = 0; i < numCoeff; i++) {
+      double dd = 1;
+      double pointCompl = 1.;
+      int exponentCompl = order;
+      for (int k = 0; k < dimSimplex; k++) {
+        dd *= nChoosek(exponentCompl, (int) exponents(i, k))
+          * pow(bezuvw[k], exponents(i, k));
+        pointCompl -= bezuvw[k];
+        exponentCompl -= (int) exponents(i, k);
+      }
+      dd *= pow(pointCompl, exponentCompl);
+
+      for (int k = dimSimplex; k < dim; k++)
+        dd *= nChoosek(order, (int) exponents(i, k))
+            * pow(bezuvw[k], exponents(i, k))
+            * pow(1. - bezuvw[k], order - exponents(i, k));
+      terms[t] += _coefficientsBez(i, t) * dd;
+    }
+  }
+
+  switch (_coefficientsBez.size2()) {
+  case 1:
+    minmaxQ[0] = terms[0];
+    minmaxQ[1] = terms[0];
+    break;
+
+  case 3:
+  {
+    double tmp = pow(terms[1], 2);
+    tmp += pow(terms[2], 2);
+    tmp = std::sqrt(tmp);
+    minmaxQ[0] = terms[0] - tmp;
+    minmaxQ[1] = terms[0] + tmp;
+  }
+    break;
+
+  case 7:
+  {
+    double tmp = pow(terms[1], 2);
+    tmp += pow(terms[2], 2);
+    tmp += pow(terms[3], 2);
+    tmp += 2 * pow(terms[4], 2);
+    tmp += 2 * pow(terms[5], 2);
+    tmp += 2 * pow(terms[6], 2);
+    tmp = std::sqrt(tmp);
+    double factor = std::sqrt(6)/3;
+    if (tmp < 1e-3*terms[0]) {
+      static int aa = 0;
+      Msg::Info("%d", ++aa);
+      minmaxQ[0] = terms[0] - factor * tmp;
+      minmaxQ[1] = terms[0] + factor * tmp;
+    }
+    else {
+      double phi;
+      {
+        fullMatrix<double> nodes(_jacobian->getNumMapNodes(),3);
+        _element->getNodesCoord(nodes);
+        fullVector<double> jac(_jacobian->getNumJacNodes());
+        _jacobian->getSignedJacobian(nodes, jac);
+
+        nodes.resize(1, 3);
+        nodes(0, 0) = uvw[0];
+        nodes(0, 1) = uvw[1];
+        nodes(0, 2) = uvw[2];
+
+        fullMatrix<double> result;
+        _jacobian->interpolate(jac, nodes, result);
+        phi = result(0, 0)*result(0, 0);
+      }
+      phi -= terms[0]*terms[0]*terms[0];
+      phi += .5*terms[0]*tmp*tmp;
+      phi /= tmp*tmp*tmp;
+      phi *= 3*std::sqrt(6);
+      phi = std::acos(phi)/3;
+      minmaxQ[0] = terms[0] + factor * tmp * std::cos(phi + 2*M_PI/3);
+      minmaxQ[1] = terms[0] + factor * tmp * std::cos(phi);
+    }
+  }
+    break;
+
+  default:
+    Msg::Error("Wrong number of functions for metric: %d",
+               _coefficientsLag.size2());
+  }
+}
+
+void MetricCoefficient::_interpolateBezierPyramid(const double *uvw, double *minmaxQ)
+{
+  int order = _jacobian->bezier->getOrder();
+
+  fullMatrix<double> exponents = JacobianBasis::generateJacMonomialsPyramid(order);
+  int numCoeff = exponents.size1();
+
+  if (_coefficientsBez.size2() == 0) _computeBezCoeff();
+
+  static int aa = 0;
+  if (++aa == 1) {
+    fullMatrix<double> val;
+    getCoefficients(val, true);
+    for (int i = 0; i < val.size2(); ++i) {
+      Msg::Info("--------- column %d", i);
+      for (int j = 0; j < val.size1(); ++j) {
+        Msg::Info("%.2e", val(j, i));
+      }
+    }
+  }
+  double terms[7];
+  for (int t = 0; t < _coefficientsBez.size2(); ++t) {
+    terms[t] = 0;
+    for (int j = 0; j < numCoeff; j++) {
+      double dd =
+          nChoosek(order, exponents(j, 2))
+          * pow(uvw[2], exponents(j, 2))
+          * pow(1. - uvw[2], order - exponents(j, 2));
+
+      double p = order + 2 - exponents(j, 2);
+      double denom = 1. - uvw[2];
+      dd *= nChoosek(p, exponents(j, 0))
+          * nChoosek(p, exponents(j, 1))
+          * pow(uvw[0] / denom, exponents(j, 0))
+          * pow(uvw[1] / denom, exponents(j, 1))
+          * pow(1. - uvw[0] / denom, p - exponents(j, 0))
+          * pow(1. - uvw[1] / denom, p - exponents(j, 1));
+      terms[t] += _coefficientsBez(j, t) * dd;
+    }
+  }
+
+  double tmp = pow(terms[1], 2);
+  tmp += pow(terms[2], 2);
+  tmp += pow(terms[3], 2);
+  tmp += 2 * pow(terms[4], 2);
+  tmp += 2 * pow(terms[5], 2);
+  tmp += 2 * pow(terms[6], 2);
+  tmp = std::sqrt(tmp);
+  double factor = std::sqrt(6)/3;
+  minmaxQ[0] = terms[0] - factor * tmp;
+  minmaxQ[1] = terms[0] + factor * tmp;
+}
+
+void MetricCoefficient::_computeBezCoeff()
+{
+  _coefficientsBez.resize(_coefficientsLag.size1(),
+                          _coefficientsLag.size2());
+  _jacobian->lag2Bez(_coefficientsLag, _coefficientsBez);
+}

Numeric/MetricBasis.h

+// Gmsh - Copyright (C) 1997-2013 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@geuz.org>.
+
+#ifndef _METRIC_BASIS_H_
+#define _METRIC_BASIS_H_
+
+#include "MElement.h"
+#include "JacobianBasis.h"
+#include "fullMatrix.h"
+
+class MetricCoefficient {
+ private:
+  MElement *_element;
+  const JacobianBasis *_jacobian;
+  fullMatrix<double> _coefficientsLag, _coefficientsBez;
+
+ public:
+  MetricCoefficient(MElement*);
+  void getCoefficients(fullMatrix<double>&, bool bezier = true);
+  void interpolate(const double *uvw, double *minmaxQ);
+  void getBound(double tol);
+
+ private:
+  void _computeBezCoeff();
+  void _interpolateBezierPyramid(const double *uvw, double *minmaxQ);
+};
+
+#endif

Numeric/bezierBasis.cpp

@@ -453,16 +453,84 @@ namespace {
   }
 }
 
+void bezierBasis::interpolate(const fullMatrix<double> &coeffs,
+                              const fullMatrix<double> &uvw,
+                              fullMatrix<double> &result,
+                              bool bezCoord) const
+{
+  result.resize(uvw.size1(), coeffs.size2());
+  int dimSimplex;
+  fullMatrix<double> bezuvw = uvw;
+  switch (type) {
+  case TYPE_HEX:
+    if (!bezCoord) {
+      bezuvw.add(1);
+      bezuvw.scale(.5);
+    }
+    dimSimplex = 0;
+    break;
+
+  case TYPE_TET:
+    dimSimplex = 3;
+    break;
+
+  case TYPE_PRI:
+    if (!bezCoord) {
+      fullMatrix<double> tmp;
+      tmp.setAsProxy(bezuvw, 3, 1);
+      tmp.add(1);
+      tmp.scale(.5);
+    }
+    dimSimplex = 2;
+    break;
+
+  default:
+  case TYPE_PYR:
+    Msg::Error("Bezier interpolation not implemented for type of element %d", type);
+    /*bezuvw[0] = .5 * (1 + uvw[0]);
+    bezuvw[1] = .5 * (1 + uvw[1]);
+    bezuvw[2] = uvw[2];
+    _interpolateBezierPyramid(uvw, minmaxQ);*/
+    return;
+  }
+
+  int numCoeff = _exponents.size1();
+
+  for (int m = 0; m < uvw.size1(); ++m) {
+    for (int n = 0; n < coeffs.size2(); ++n) result(m, n) = 0;
+    for (int i = 0; i < numCoeff; i++) {
+      double dd = 1;
+      double pointCompl = 1.;
+      int exponentCompl = order;
+      for (int k = 0; k < dimSimplex; k++) {
+        dd *= nChoosek(exponentCompl, (int) _exponents(i, k))
+          * pow(bezuvw(m, k), _exponents(i, k));
+        pointCompl -= bezuvw(m, k);
+        exponentCompl -= (int) _exponents(i, k);
+      }
+      dd *= pow_int(pointCompl, exponentCompl);
+
+      for (int k = dimSimplex; k < dim; k++) {
+        dd *= nChoosek(order, (int) _exponents(i, k))
+            * pow_int(bezuvw(m, k), _exponents(i, k))
+            * pow_int(1. - bezuvw(m, k), order - _exponents(i, k));
+      }
+      for (int n = 0; n < coeffs.size2(); ++n)
+        result(m, n) += coeffs(i, n) * dd;
+    }
+  }
+}
+
 void bezierBasis::_construct(int parentType, int p)
 {
   order = p;
-  fullMatrix<double> exponents;
+  type = parentType;
   std::vector< fullMatrix<double> > subPoints;
 
   if (parentType == TYPE_PYR) {
     dim = 3;
     numLagCoeff = 8;
-    exponents = JacobianBasis::generateJacMonomialsPyramid(order);
+    _exponents = JacobianBasis::generateJacMonomialsPyramid(order);
     if (order < 2) {
       subPoints = generateSubPointsPyr(order);
       numDivisions = static_cast<int>(subPoints.size());
@@ -473,13 +541,13 @@ void bezierBasis::_construct(int parentType, int p)
       if (++a == 1) Msg::Warning("Subdivision not available for pyramid p>1");
     }
 
-    fullMatrix<double> bezierPoints = exponents;
+    fullMatrix<double> bezierPoints = _exponents;
     bezierPoints.scale(1. / (order + 2));
 
-    matrixBez2Lag = generateBez2LagMatrixPyramid(exponents, bezierPoints, order);
+    matrixBez2Lag = generateBez2LagMatrixPyramid(_exponents, bezierPoints, order);
     matrixBez2Lag.invert(matrixLag2Bez);
     if (order < 2)
-      subDivisor = generateSubDivisorPyramid(exponents, subPoints, matrixLag2Bez, order);
+      subDivisor = generateSubDivisorPyramid(_exponents, subPoints, matrixLag2Bez, order);
     return;
   }
 
@@ -489,14 +557,14 @@ void bezierBasis::_construct(int parentType, int p)
       dim = 0;
       numLagCoeff = 1;
       dimSimplex = 0;
-      exponents = gmshGenerateMonomialsLine(0);
+      _exponents = gmshGenerateMonomialsLine(0);
       subPoints.push_back(gmshGeneratePointsLine(0));
       break;
     case TYPE_LIN : {
       dim = 1;
       numLagCoeff = 2;
       dimSimplex = 0;
-      exponents = gmshGenerateMonomialsLine(order);
+      _exponents = gmshGenerateMonomialsLine(order);
       subPoints = generateSubPointsLine(order);
       break;
     }
@@ -504,7 +572,7 @@ void bezierBasis::_construct(int parentType, int p)
       dim = 2;
       numLagCoeff = 3;
       dimSimplex = 2;
-      exponents = gmshGenerateMonomialsTriangle(order);
+      _exponents = gmshGenerateMonomialsTriangle(order);
       subPoints = generateSubPointsTriangle(order);
       break;
     }
@@ -512,7 +580,7 @@ void bezierBasis::_construct(int parentType, int p)
       dim = 2;
       numLagCoeff = 4;
       dimSimplex = 0;
-      exponents = gmshGenerateMonomialsQuadrangle(order);
+      _exponents = gmshGenerateMonomialsQuadrangle(order);
       subPoints = generateSubPointsQuad(order);
       break;
     }
@@ -520,7 +588,7 @@ void bezierBasis::_construct(int parentType, int p)
       dim = 3;
       numLagCoeff = 4;
       dimSimplex = 3;
-      exponents = gmshGenerateMonomialsTetrahedron(order);
+      _exponents = gmshGenerateMonomialsTetrahedron(order);
       subPoints = generateSubPointsTetrahedron(order);
       break;
     }
@@ -528,7 +596,7 @@ void bezierBasis::_construct(int parentType, int p)
       dim = 3;
       numLagCoeff = 6;
       dimSimplex = 2;
-      exponents = gmshGenerateMonomialsPrism(order);
+      _exponents = gmshGenerateMonomialsPrism(order);
       subPoints = generateSubPointsPrism(order);
       break;
     }
@@ -536,7 +604,7 @@ void bezierBasis::_construct(int parentType, int p)
       dim = 3;
       numLagCoeff = 8;
       dimSimplex = 0;
-      exponents = gmshGenerateMonomialsHexahedron(order);
+      _exponents = gmshGenerateMonomialsHexahedron(order);
       subPoints = generateSubPointsHex(order);
       break;
     }
@@ -547,17 +615,17 @@ void bezierBasis::_construct(int parentType, int p)
       order = 0;
       numLagCoeff = 4;
       dimSimplex = 3;
-      exponents = gmshGenerateMonomialsTetrahedron(order);
+      _exponents = gmshGenerateMonomialsTetrahedron(order);
       subPoints = generateSubPointsTetrahedron(order);
       break;
     }
   }
   numDivisions = static_cast<int>(subPoints.size());
 
-  fullMatrix<double> bezierPoints = exponents;
+  fullMatrix<double> bezierPoints = _exponents;
   bezierPoints.scale(1./order);
 
-  matrixBez2Lag = generateBez2LagMatrix(exponents, bezierPoints, order, dimSimplex);
+  matrixBez2Lag = generateBez2LagMatrix(_exponents, bezierPoints, order, dimSimplex);
   matrixBez2Lag.invert(matrixLag2Bez);
-  subDivisor = generateSubDivisor(exponents, subPoints, matrixLag2Bez, order, dimSimplex);
+  subDivisor = generateSubDivisor(_exponents, subPoints, matrixLag2Bez, order, dimSimplex);
 }

Numeric/bezierBasis.h

@@ -17,8 +17,9 @@ class bezierBasis {
  private :
   // the 'numLagCoeff' first exponents are related to 'Lagrangian' values
   int numLagCoeff;
-  int dim, order;
+  int dim, type, order;
   int numDivisions;
+  fullMatrix<double> _exponents;
 
  public :
   fullMatrix<double> matrixLag2Bez;
@@ -39,6 +40,14 @@ class bezierBasis {
   inline int getNumLagCoeff() const {return numLagCoeff;}
   inline int getNumDivision() const {return numDivisions;}
 
+  // Interpolation of n functions on N points :
+  // coeffs(numCoeff, n) and uvw(N, dim)
+  // => result(N, n)
+  void interpolate(const fullMatrix<double> &coeffs,
+                   const fullMatrix<double> &uvw,
+                   fullMatrix<double> &result,
+                   bool bezCoord = false) const;
+
  private :
   void _construct(int parendtType, int order);
 };