diff --git a/Numeric/BasisFactory.cpp b/Numeric/BasisFactory.cpp
index 129011088747e34b91fab94e7efcc5c1a119990a..96f51530117c1f88b6997d54ba73923cf34b87b9 100644
--- a/Numeric/BasisFactory.cpp
+++ b/Numeric/BasisFactory.cpp
@@ -14,6 +14,7 @@
std::map<int, nodalBasis*> BasisFactory::fs;
std::map<int, JacobianBasis*> BasisFactory::js;
BasisFactory::Cont_bezierBasis BasisFactory::bs;
+BasisFactory::Cont_gradBasis BasisFactory::gs;
const nodalBasis* BasisFactory::getNodalBasis(int tag)
{
@@ -69,14 +70,27 @@ const JacobianBasis* BasisFactory::getJacobianBasis(int tag)
return J;
}
+const GradientBasis* BasisFactory::getGradientBasis(int tag, int order)
+{
+ std::pair<int, int> key(tag, order);
+ Cont_gradBasis::const_iterator it = gs.find(key);
+ if (it != gs.end())
+ return it->second;
+
+ GradientBasis* G = new GradientBasis(tag, order);
+ gs.insert(std::make_pair(key, G));
+ return G;
+}
+
const bezierBasis* BasisFactory::getBezierBasis(int parentType, int order)
{
- Cont_bezierBasis::iterator it = bs.find(std::make_pair(parentType, order));
+ std::pair<int, int> key(parentType, order);
+ Cont_bezierBasis::iterator it = bs.find(key);
if (it != bs.end())
return it->second;
bezierBasis* B = new bezierBasis(parentType, order);
- bs.insert(std::make_pair(std::make_pair(parentType, order), B));
+ bs.insert(std::make_pair(key, B));
return B;
}
diff --git a/Numeric/BasisFactory.h b/Numeric/BasisFactory.h
index d7b8e51fcfd3d5c771fcb96cf7e860d836862995..171f740c340f36d1001cad3b4217561868d59ad3 100644
--- a/Numeric/BasisFactory.h
+++ b/Numeric/BasisFactory.h
@@ -14,10 +14,12 @@
class BasisFactory
{
typedef std::map<std::pair<int, int>, bezierBasis*> Cont_bezierBasis;
+ typedef std::map<std::pair<int, int>, GradientBasis*> Cont_gradBasis;
private:
static std::map<int, nodalBasis*> fs;
static std::map<int, JacobianBasis*> js;
+ static Cont_gradBasis gs;
static Cont_bezierBasis bs;
// store bezier bases by parentType and order (no serendipity..)
@@ -25,6 +27,7 @@ class BasisFactory
// Caution: the returned pointer can be NULL
static const nodalBasis* getNodalBasis(int tag);
static const JacobianBasis* getJacobianBasis(int tag);
+ static const GradientBasis* getGradientBasis(int tag, int order);
static const bezierBasis* getBezierBasis(int parentTag, int order);
static inline const bezierBasis* getBezierBasis(int tag) {
return getBezierBasis(ElementType::ParentTypeFromTag(tag),
diff --git a/Numeric/JacobianBasis.cpp b/Numeric/JacobianBasis.cpp
index 40a18eac4f094244fabcbc4ef391f32c79797050..f39583a4e532ee971dcc97880682dfab16e8aa02 100644
--- a/Numeric/JacobianBasis.cpp
+++ b/Numeric/JacobianBasis.cpp
@@ -23,6 +23,73 @@ namespace {
}
}
+GradientBasis::GradientBasis(int tag, int order)
+{
+ _bezier = BasisFactory::getBezierBasis(ElementType::ParentTypeFromTag(tag), order);
+ const int type = ElementType::ParentTypeFromTag(tag);
+
+ fullMatrix<double> samplingPoints;
+
+ switch (type) {
+ case TYPE_PNT :
+ samplingPoints = gmshGeneratePointsLine(0);
+ break;
+ case TYPE_LIN :
+ samplingPoints = gmshGeneratePointsLine(order);
+ break;
+ case TYPE_TRI :
+ samplingPoints = gmshGeneratePointsTriangle(order,false);
+ break;
+ case TYPE_QUA :
+ samplingPoints = gmshGeneratePointsQuadrangle(order,false);
+ break;
+ case TYPE_TET :
+ samplingPoints = gmshGeneratePointsTetrahedron(order,false);
+ break;
+ case TYPE_PRI :
+ samplingPoints = gmshGeneratePointsPrism(order,false);
+ break;
+ case TYPE_HEX :
+ samplingPoints = gmshGeneratePointsHexahedron(order,false);
+ break;
+ case TYPE_PYR :
+ Msg::Error("not sure if pyramidal metric space = jacobian space");
+ //samplingPoints = JacobianBasis::generateJacPointsPyramid(order);
+ break;
+ default :
+ Msg::Error("Unknown Jacobian function space for element type %d", type);
+ return;
+ }
+ const int numSampPnts = samplingPoints.size1();
+
+ // Store shape function gradients of mapping at Jacobian nodes
+ fullMatrix<double> allDPsi;
+ const nodalBasis *mapBasis = BasisFactory::getNodalBasis(tag);
+ mapBasis->df(samplingPoints, allDPsi);
+ const int numMapNodes = allDPsi.size1();
+
+ gradShapeMatX.resize(numSampPnts, numMapNodes);
+ gradShapeMatY.resize(numSampPnts, numMapNodes);
+ gradShapeMatZ.resize(numSampPnts, numMapNodes);
+ for (int i=0; i<numSampPnts; i++) {
+ for (int j=0; j<numMapNodes; j++) {
+ gradShapeMatX(i, j) = allDPsi(j, 3*i);
+ gradShapeMatY(i, j) = allDPsi(j, 3*i+1);
+ gradShapeMatZ(i, j) = allDPsi(j, 3*i+2);
+ }
+ }
+}
+
+void GradientBasis::getGradientsFromNodes(const fullMatrix<double> &nodes,
+ fullMatrix<double> *dxyzdX,
+ fullMatrix<double> *dxyzdY,
+ fullMatrix<double> *dxyzdZ) const
+{
+ if (dxyzdX) gradShapeMatX.mult(nodes, *dxyzdX);
+ if (dxyzdY) gradShapeMatY.mult(nodes, *dxyzdY);
+ if (dxyzdZ) gradShapeMatZ.mult(nodes, *dxyzdZ);
+}
+
JacobianBasis::JacobianBasis(int tag)
{
const int parentType = ElementType::ParentTypeFromTag(tag);
@@ -58,7 +125,7 @@ JacobianBasis::JacobianBasis(int tag)
lagPoints = generateJacPointsPyramid(jacobianOrder);
break;
default :
- Msg::Error("Unknown Jacobian function space for element type %d", tag);
+ Msg::Error("Unknown Jacobian function space for element type %d", parentType);
return;
}
numJacNodes = lagPoints.size1();
@@ -67,22 +134,7 @@ JacobianBasis::JacobianBasis(int tag)
bezier = BasisFactory::getBezierBasis(parentType, jacobianOrder);
// Store shape function gradients of mapping at Jacobian nodes
- fullMatrix<double> allDPsi;
- const nodalBasis *mapBasis = BasisFactory::getNodalBasis(tag);
-
- mapBasis->df(lagPoints, allDPsi);
- numMapNodes = allDPsi.size1();
-
- gradShapeMatX.resize(numJacNodes, numMapNodes);
- gradShapeMatY.resize(numJacNodes, numMapNodes);
- gradShapeMatZ.resize(numJacNodes, numMapNodes);
- for (int i=0; i<numJacNodes; i++) {
- for (int j=0; j<numMapNodes; j++) {
- gradShapeMatX(i, j) = allDPsi(j, 3*i);
- gradShapeMatY(i, j) = allDPsi(j, 3*i+1);
- gradShapeMatZ(i, j) = allDPsi(j, 3*i+2);
- }
- }
+ _gradBasis = new GradientBasis(tag, jacobianOrder);
// Compute matrix for lifting from primary Jacobian basis to Jacobian basis
int primJacType = ElementType::getTag(parentType, primJacobianOrder, false);
@@ -133,7 +185,9 @@ JacobianBasis::JacobianBasis(int tag)
lagPointsFast(numPrimMapNodes,2) = barycenter[2];
fullMatrix<double> allDPsiFast;
+ const nodalBasis *mapBasis = BasisFactory::getNodalBasis(tag);
mapBasis->df(lagPointsFast, allDPsiFast);
+ numMapNodes = mapBasis->getNumShapeFunctions();
gradShapeMatXFast.resize(numJacNodesFast, numMapNodes);
gradShapeMatYFast.resize(numJacNodesFast, numMapNodes);
@@ -572,8 +626,8 @@ void JacobianBasis::getMetricMinAndGradients(const fullMatrix<double> &nodesXYZ,
for (int l = 0; l < numJacNodes; l++) {
double jac[2][2] = {{0., 0.}, {0., 0.}};
for (int i = 0; i < numMapNodes; i++) {
- const double &dPhidX = gradShapeMatX(l,i);
- const double &dPhidY = gradShapeMatY(l,i);
+ const double &dPhidX = _gradBasis->getgradShapeMatX(l,i);
+ const double &dPhidY = _gradBasis->getgradShapeMatY(l,i);
const double dpsidx = dPhidX * invJaci[0][0] + dPhidY * invJaci[1][0];
const double dpsidy = dPhidX * invJaci[0][1] + dPhidY * invJaci[1][1];
jac[0][0] += nodesXYZ(i,0) * dpsidx;
@@ -596,8 +650,8 @@ void JacobianBasis::getMetricMinAndGradients(const fullMatrix<double> &nodesXYZ,
const double aetaxi = ayx * invJaci[0][0] + ayy * invJaci[0][1];
const double axieta = axx * invJaci[1][0] + axy * invJaci[1][1];
for (int i = 0; i < numMapNodes; i++) {
- const double &dPhidX = gradShapeMatX(l,i);
- const double &dPhidY = gradShapeMatY(l,i);
+ const double &dPhidX = _gradBasis->gradShapeMatX(l,i);
+ const double &dPhidY = _gradBasis->gradShapeMatY(l,i);
gradLambdaJ(l, i + 0 * numMapNodes) = axixi * dPhidX + axieta * dPhidY;
gradLambdaJ(l, i + 1 * numMapNodes) = aetaxi * dPhidX + aetaeta * dPhidY;
}
@@ -759,12 +813,3 @@ int JacobianBasis::jacobianOrder(int parentType, int order)
}
-void JacobianBasis::getGradientsFromNodes(const fullMatrix<double> &nodes,
- fullMatrix<double> *dxyzdX,
- fullMatrix<double> *dxyzdY,
- fullMatrix<double> *dxyzdZ) const
-{
- if (dxyzdX) gradShapeMatX.mult(nodes, *dxyzdX);
- if (dxyzdY) gradShapeMatY.mult(nodes, *dxyzdY);
- if (dxyzdZ) gradShapeMatZ.mult(nodes, *dxyzdZ);
-}
diff --git a/Numeric/JacobianBasis.h b/Numeric/JacobianBasis.h
index 2c71d649bdc8dec293988bc80b94b7f2e4a6a857..312808ab890eb5f65b6807f2fe7a1ba23b1974bd 100644
--- a/Numeric/JacobianBasis.h
+++ b/Numeric/JacobianBasis.h
@@ -12,9 +12,40 @@
#include "fullMatrix.h"
+class GradientBasis {
+ public:
+ fullMatrix<double> gradShapeMatX, gradShapeMatY, gradShapeMatZ;
+
+ private:
+ const bezierBasis *_bezier;
+
+ public :
+
+ GradientBasis(int tag, int order);
+
+ int getNumSamplingPoints() const {return gradShapeMatX.size1();}
+ int getNumMapNodes() const {return gradShapeMatX.size2();}
+ const bezierBasis* getBezierBasis() const {return _bezier;}
+
+ void getGradientsFromNodes(const fullMatrix<double> &nodes,
+ fullMatrix<double> *dxyzdX,
+ fullMatrix<double> *dxyzdY,
+ fullMatrix<double> *dxyzdZ) const;
+ inline double getgradShapeMatX(int i, int j) const {return gradShapeMatX(i, j);}
+ inline double getgradShapeMatY(int i, int j) const {return gradShapeMatY(i, j);}
+ inline double getgradShapeMatZ(int i, int j) const {return gradShapeMatZ(i, j);}
+ inline void lag2Bez(const fullVector<double> &jac, fullVector<double> &bez) const {
+ _bezier->matrixLag2Bez.mult(jac,bez);
+ }
+ inline void lag2Bez(const fullMatrix<double> &jac, fullMatrix<double> &bez) const {
+ _bezier->matrixLag2Bez.mult(jac,bez);
+ }
+};
+
+
class JacobianBasis {
private:
- fullMatrix<double> gradShapeMatX, gradShapeMatY, gradShapeMatZ;
+ GradientBasis *_gradBasis;
fullMatrix<double> gradShapeMatXFast, gradShapeMatYFast, gradShapeMatZFast;
fullVector<double> primGradShapeBarycenterX, primGradShapeBarycenterY, primGradShapeBarycenterZ;
fullMatrix<double> matrixPrimJac2Jac; // Lifts Lagrange basis of primary Jac. to Lagrange basis of Jac.
@@ -23,21 +54,6 @@ class JacobianBasis {
int numMapNodes, numPrimMapNodes;
int numJacNodesFast;
- void getSignedJacobianGeneral(int nJacNodes, const fullMatrix<double> &gSMatX,
- const fullMatrix<double> &gSMatY, const fullMatrix<double> &gSMatZ,
- const fullMatrix<double> &nodesXYZ, fullVector<double> &jacobian) const;
- void getSignedJacobianGeneral(int nJacNodes, const fullMatrix<double> &gSMatX,
- const fullMatrix<double> &gSMatY, const fullMatrix<double> &gSMatZ,
- const fullMatrix<double> &nodesX, const fullMatrix<double> &nodesY,
- const fullMatrix<double> &nodesZ, fullMatrix<double> &jacobian) const;
- void getScaledJacobianGeneral(int nJacNodes, const fullMatrix<double> &gSMatX,
- const fullMatrix<double> &gSMatY, const fullMatrix<double> &gSMatZ,
- const fullMatrix<double> &nodesXYZ, fullVector<double> &jacobian) const;
- void getSignedJacAndGradientsGeneral(int nJacNodes, const fullMatrix<double> &gSMatX,
- const fullMatrix<double> &gSMatY, const fullMatrix<double> &gSMatZ,
- const fullMatrix<double> &nodesXYZ, const fullMatrix<double> &normals,
- fullMatrix<double> &JDJ) const;
-
public :
const bezierBasis *bezier;
@@ -59,7 +75,8 @@ class JacobianBasis {
double getPrimJac3D(const fullMatrix<double> &nodesXYZ) const;
inline void getSignedJacAndGradients(const fullMatrix<double> &nodesXYZ,
const fullMatrix<double> &normals, fullMatrix<double> &JDJ) const {
- getSignedJacAndGradientsGeneral(numJacNodes,gradShapeMatX,gradShapeMatY,gradShapeMatZ,nodesXYZ,normals,JDJ);
+ getSignedJacAndGradientsGeneral(numJacNodes, _gradBasis->gradShapeMatX,
+ _gradBasis->gradShapeMatY, _gradBasis->gradShapeMatZ, nodesXYZ, normals, JDJ);
}
inline void getSignedJacAndGradientsFast(const fullMatrix<double> &nodesXYZ,
const fullMatrix<double> &normals, fullMatrix<double> &JDJ) const {
@@ -70,15 +87,17 @@ class JacobianBasis {
const fullMatrix<double> &nodesXYZStraight,
fullVector<double> &lambdaJ , fullMatrix<double> &gradLambdaJ) const;
inline void getSignedJacobian(const fullMatrix<double> &nodesXYZ, fullVector<double> &jacobian) const {
- getSignedJacobianGeneral(numJacNodes,gradShapeMatX,gradShapeMatY,gradShapeMatZ,nodesXYZ,jacobian);
+ getSignedJacobianGeneral(numJacNodes, _gradBasis->gradShapeMatX,
+ _gradBasis->gradShapeMatY, _gradBasis->gradShapeMatZ,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);
+ getSignedJacobianGeneral(numJacNodes, _gradBasis->gradShapeMatX,
+ _gradBasis->gradShapeMatY, _gradBasis->gradShapeMatZ,nodesX,nodesY,nodesZ,jacobian);
}
inline void getScaledJacobian(const fullMatrix<double> &nodesXYZ, fullVector<double> &jacobian) const {
- getScaledJacobianGeneral(numJacNodes,gradShapeMatX,gradShapeMatY,gradShapeMatZ,nodesXYZ,jacobian);
+ getScaledJacobianGeneral(numJacNodes, _gradBasis->gradShapeMatX,
+ _gradBasis->gradShapeMatY, _gradBasis->gradShapeMatZ,nodesXYZ,jacobian);
}
inline void getSignedJacobianFast(const fullMatrix<double> &nodesXYZ, fullVector<double> &jacobian) const {
getSignedJacobianGeneral(numJacNodesFast,gradShapeMatXFast,gradShapeMatYFast,gradShapeMatZFast,nodesXYZ,jacobian);
@@ -109,16 +128,27 @@ class JacobianBasis {
const fullMatrix<double> &uvw,
fullMatrix<double> &result) const;
- // Jacobian basis order and pyramidal basis
- void getGradientsFromNodes(const fullMatrix<double> &nodes,
- fullMatrix<double> *dxyzdX,
- fullMatrix<double> *dxyzdY,
- fullMatrix<double> *dxyzdZ) const;
-
//
static int jacobianOrder(int parentType, int order);
static fullMatrix<double> generateJacMonomialsPyramid(int order);
static fullMatrix<double> generateJacPointsPyramid(int order);
+
+
+ private :
+ void getSignedJacobianGeneral(int nJacNodes, const fullMatrix<double> &gSMatX,
+ const fullMatrix<double> &gSMatY, const fullMatrix<double> &gSMatZ,
+ const fullMatrix<double> &nodesXYZ, fullVector<double> &jacobian) const;
+ void getSignedJacobianGeneral(int nJacNodes, const fullMatrix<double> &gSMatX,
+ const fullMatrix<double> &gSMatY, const fullMatrix<double> &gSMatZ,
+ const fullMatrix<double> &nodesX, const fullMatrix<double> &nodesY,
+ const fullMatrix<double> &nodesZ, fullMatrix<double> &jacobian) const;
+ void getScaledJacobianGeneral(int nJacNodes, const fullMatrix<double> &gSMatX,
+ const fullMatrix<double> &gSMatY, const fullMatrix<double> &gSMatZ,
+ const fullMatrix<double> &nodesXYZ, fullVector<double> &jacobian) const;
+ void getSignedJacAndGradientsGeneral(int nJacNodes, const fullMatrix<double> &gSMatX,
+ const fullMatrix<double> &gSMatY, const fullMatrix<double> &gSMatZ,
+ const fullMatrix<double> &nodesXYZ, const fullMatrix<double> &normals,
+ fullMatrix<double> &JDJ) const;
};
#endif
diff --git a/Numeric/MetricBasis.cpp b/Numeric/MetricBasis.cpp
index f1b2292d8a06f296ab9b342158f60e4286e1a48c..865f31aee3e7d71290fd858761cd7f042c6228d4 100644
--- a/Numeric/MetricBasis.cpp
+++ b/Numeric/MetricBasis.cpp
@@ -11,10 +11,13 @@
MetricCoefficient::MetricCoefficient(MElement *el) : _element(el)
{
- _jacobian = BasisFactory::getJacobianBasis(el->getTypeForMSH());
+ const int tag = el->getTypeForMSH();
+ const int metricOrder = MetricCoefficient::metricOrder(tag);
+ _jacobian = BasisFactory::getJacobianBasis(tag);
+ _gradients = BasisFactory::getGradientBasis(tag, metricOrder);
- int nJac = _jacobian->getNumJacNodes();
- int nMapping = _jacobian->getNumMapNodes();
+ int nSampPnts = _gradients->getNumSamplingPoints();
+ int nMapping = _gradients->getNumMapNodes();
fullMatrix<double> nodesXYZ(nMapping, 3);
el->getNodesCoord(nodesXYZ);
@@ -36,11 +39,11 @@ MetricCoefficient::MetricCoefficient(MElement *el) : _element(el)
case 3 :
{
- fullMatrix<double> dxyzdX(nJac,3), dxyzdY(nJac,3), dxyzdZ(nJac,3);
- _jacobian->getGradientsFromNodes(nodesXYZ, &dxyzdX, &dxyzdY, &dxyzdZ);
+ fullMatrix<double> dxyzdX(nSampPnts,3), dxyzdY(nSampPnts,3), dxyzdZ(nSampPnts,3);
+ _gradients->getGradientsFromNodes(nodesXYZ, &dxyzdX, &dxyzdY, &dxyzdZ);
- _coefficientsLag.resize(nJac, 7);
- for (int i = 0; i < nJac; i++) {
+ _coefficientsLag.resize(nSampPnts, 7);
+ for (int i = 0; i < nSampPnts; 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);
@@ -140,7 +143,7 @@ void MetricCoefficient::interpolate(const double *uvw, double *minmaxQ)
return;
}
- int order = _jacobian->bezier->getOrder();
+ int order = _gradients->getBezierBasis()->getOrder();
int dimSimplex;
fullMatrix<double> exponents;
@@ -241,7 +244,7 @@ void MetricCoefficient::interpolate(const double *uvw, double *minmaxQ)
else {
double phi;
{
- fullMatrix<double> nodes(_jacobian->getNumMapNodes(),3);
+ fullMatrix<double> nodes(_element->getNumShapeFunctions(),3);
_element->getNodesCoord(nodes);
fullVector<double> jac(_jacobian->getNumJacNodes());
_jacobian->getSignedJacobian(nodes, jac);
@@ -259,6 +262,8 @@ void MetricCoefficient::interpolate(const double *uvw, double *minmaxQ)
phi += .5*terms[0]*tmp*tmp;
phi /= tmp*tmp*tmp;
phi *= 3*std::sqrt(6);
+ if (phi > 1) phi = 1;
+ if (phi < -1) phi = -1;
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);
@@ -272,6 +277,27 @@ void MetricCoefficient::interpolate(const double *uvw, double *minmaxQ)
}
}
+
+int MetricCoefficient::metricOrder(int tag)
+{
+ const int parentType = ElementType::ParentTypeFromTag(tag);
+ const int order = ElementType::OrderFromTag(tag);
+
+ switch (parentType) {
+ case TYPE_PNT : return 0;
+ case TYPE_LIN : return order;
+ case TYPE_TRI :
+ case TYPE_TET : return 2*order-2;
+ case TYPE_QUA :
+ case TYPE_PRI :
+ case TYPE_HEX :
+ case TYPE_PYR : return 2*order;
+ default :
+ Msg::Error("Unknown element type %d, return order 0", parentType);
+ return 0;
+ }
+}
+
void MetricCoefficient::_interpolateBezierPyramid(const double *uvw, double *minmaxQ)
{
int order = _jacobian->bezier->getOrder();
@@ -329,5 +355,5 @@ void MetricCoefficient::_computeBezCoeff()
{
_coefficientsBez.resize(_coefficientsLag.size1(),
_coefficientsLag.size2());
- _jacobian->lag2Bez(_coefficientsLag, _coefficientsBez);
+ _gradients->lag2Bez(_coefficientsLag, _coefficientsBez);
}
diff --git a/Numeric/MetricBasis.h b/Numeric/MetricBasis.h
index 5c2ee6da791719262d2163a184dbc6084e485637..d1cc8a98fc29d3f3fcdc825e3047ebf6feb8329d 100644
--- a/Numeric/MetricBasis.h
+++ b/Numeric/MetricBasis.h
@@ -14,13 +14,17 @@ class MetricCoefficient {
private:
MElement *_element;
const JacobianBasis *_jacobian;
+ const GradientBasis *_gradients;
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);
+ void getBoundRmin(double tol);
+
+ static int metricOrder(int tag);
private:
void _computeBezCoeff();