diff --git a/Numeric/CMakeLists.txt b/Numeric/CMakeLists.txt
index 8abf9bc3f3749692942d81cd77b6232bae92901b..d7a0f85162762cad918498d8e384a620a25851f9 100644
--- a/Numeric/CMakeLists.txt
+++ b/Numeric/CMakeLists.txt
@@ -16,6 +16,7 @@ set(SRC
legendrePolynomials.cpp
bezierBasis.cpp
JacobianBasis.cpp
+ MetricBasis.cpp
pointsGenerators.cpp
ElementType.cpp
GaussIntegration.cpp
diff --git a/Numeric/JacobianBasis.cpp b/Numeric/JacobianBasis.cpp
index 20af53c78f1f85203916b7a58c3c8eb028e2d4a6..40a18eac4f094244fabcbc4ef391f32c79797050 100644
--- a/Numeric/JacobianBasis.cpp
+++ b/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;
diff --git a/Numeric/JacobianBasis.h b/Numeric/JacobianBasis.h
index 8a08f18e69ef38e11b2d73acaca3d1b8cfff4b95..2c71d649bdc8dec293988bc80b94b7f2e4a6a857 100644
--- a/Numeric/JacobianBasis.h
+++ b/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,
diff --git a/Numeric/MetricBasis.cpp b/Numeric/MetricBasis.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..f1b2292d8a06f296ab9b342158f60e4286e1a48c
--- /dev/null
+++ b/Numeric/MetricBasis.cpp
@@ -0,0 +1,333 @@
+// 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);
+}
diff --git a/Numeric/MetricBasis.h b/Numeric/MetricBasis.h
new file mode 100644
index 0000000000000000000000000000000000000000..5c2ee6da791719262d2163a184dbc6084e485637
--- /dev/null
+++ b/Numeric/MetricBasis.h
@@ -0,0 +1,30 @@
+// 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
diff --git a/Numeric/bezierBasis.cpp b/Numeric/bezierBasis.cpp
index 054fdaa958ddff38f2f7265284e44e1dd11ae5fe..3bc3f9c50c6aa2196deb00ac6fda6f9f475299ee 100644
--- a/Numeric/bezierBasis.cpp
+++ b/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);
}
diff --git a/Numeric/bezierBasis.h b/Numeric/bezierBasis.h
index 3ecc1784801d08453b443e5f76b1a5d7e353f153..3e236f1aae2ca2a0724bc8d35bec42686b533d29 100644
--- a/Numeric/bezierBasis.h
+++ b/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);
};