diff --git a/Numeric/MetricBasis.cpp b/Numeric/MetricBasis.cpp
index c5701a234883050525fa47de553eabad304e23ae..189a5d29535484450f424bc8b42081f5f4d79ef8 100644
--- a/Numeric/MetricBasis.cpp
+++ b/Numeric/MetricBasis.cpp
@@ -88,23 +88,6 @@ double MetricBasis::minRCorner(MElement *el)
// Nodes
fullMatrix<double> nodes(nMapping, 3);
el->getNodesCoord(nodes);
- fullMatrix<double> nodesMetric;
- switch (el->getDim()) {
- case 1:
- return -1;
-
- case 2:
- if (_paramOnPlane(nodes, nodesMetric)) return -1;
- break;
-
- case 3:
- nodesMetric.setAsProxy(nodes);
- break;
-
- default:
- Msg::Error("no metric for element of dimension %d (tag=%d)",
- el->getDim(), el->getTypeForMSH());
- }
// Jacobian coefficients
fullVector<double> jacLag(jacobian->getNumJacNodes());
@@ -112,7 +95,7 @@ double MetricBasis::minRCorner(MElement *el)
// Metric coefficients
fullMatrix<double> metCoeffLag;
- _fillCoeff(gradients, nodesMetric, metCoeffLag);
+ _fillCoeff(el->getDim(), gradients, nodes, metCoeffLag);
// Compute min_corner(R)
return _computeMinlagR(jacLag, metCoeffLag, nSampPnts);
@@ -163,16 +146,18 @@ double MetricBasis::getMinSampledR(MElement *el, int deg) const
double uvw[3];
double minmaxQ[2];
+ int dim = el->getDim();
uvw[0] = samplingPoints(0, 0);
uvw[1] = samplingPoints(0, 1);
- uvw[2] = samplingPoints(0, 2);
-
+ if (dim == 3) uvw[2] = samplingPoints(0, 2);
+ else uvw[2] = 0;
interpolate(el, md, uvw, minmaxQ);
+
double min, max = min = std::sqrt(minmaxQ[0]/minmaxQ[1]);
for (int i = 1; i < samplingPoints.size1(); ++i) {
uvw[0] = samplingPoints(i, 0);
uvw[1] = samplingPoints(i, 1);
- uvw[2] = samplingPoints(i, 2);
+ if (dim == 3) uvw[2] = samplingPoints(i, 2);
interpolate(el, md, uvw, minmaxQ);
double tmp = std::sqrt(minmaxQ[0]/minmaxQ[1]);
min = std::min(min, tmp);
@@ -189,23 +174,6 @@ double MetricBasis::getBoundMinR(MElement *el)
// Nodes
fullMatrix<double> nodes(nMapping, 3);
el->getNodesCoord(nodes);
- fullMatrix<double> nodesMetric;
- switch (el->getDim()) {
- case 1:
- return -1;
-
- case 2:
- if (_paramOnPlane(nodes, nodesMetric)) return -1;
- break;
-
- case 3:
- nodesMetric.setAsProxy(nodes);
- break;
-
- default:
- Msg::Error("no metric for element of dimension %d (tag=%d)",
- el->getDim(), el->getTypeForMSH());
- }
// Jacobian coefficients
fullVector<double> jacLag(_jacobian->getNumJacNodes());
@@ -215,7 +183,7 @@ double MetricBasis::getBoundMinR(MElement *el)
// Metric coefficients
fullMatrix<double> metCoeffLag;
- _fillCoeff(_gradients, nodesMetric, metCoeffLag);
+ _fillCoeff(el->getDim(), _gradients, nodes, metCoeffLag);
fullMatrix<double> *metCoeff;
metCoeff = new fullMatrix<double>(nSampPnts, metCoeffLag.size2());
_bezier->matrixLag2Bez.mult(metCoeffLag, *metCoeff);
@@ -453,6 +421,14 @@ void MetricBasis::interpolate(const MElement *el, const MetricData *md, const do
dimSimplex = 2;
exponents = gmshGenerateMonomialsPrism(order);
break;
+
+ case TYPE_TRI:
+ bezuvw[0] = uvw[0];
+ bezuvw[1] = uvw[1];
+ bezuvw[2] = uvw[2];
+ dimSimplex = 2;
+ exponents = gmshGenerateMonomialsTriangle(order);
+ break;
}
int numCoeff = exponents.size1();
@@ -823,27 +799,6 @@ void MetricBasis::_getMetricData(MElement *el, MetricData *&md) const
// Nodes
fullMatrix<double> nodes(nMapping, 3);
el->getNodesCoord(nodes);
- fullMatrix<double> nodesMetric;
- switch (el->getDim()) {
- case 1:
- md = NULL;
- return;
-
- case 2:
- if (_paramOnPlane(nodes, nodesMetric)) {
- md = NULL;
- return;
- }
- break;
-
- case 3:
- nodesMetric.setAsProxy(nodes);
- break;
-
- default:
- Msg::Error("no metric for element of dimension %d (tag=%d)",
- el->getDim(), el->getTypeForMSH());
- }
// Jacobian coefficients
fullVector<double> jacLag(_jacobian->getNumJacNodes());
@@ -853,7 +808,7 @@ void MetricBasis::_getMetricData(MElement *el, MetricData *&md) const
// Metric coefficients
fullMatrix<double> metCoeffLag;
- _fillCoeff(_gradients, nodesMetric, metCoeffLag);
+ _fillCoeff(el->getDim(), _gradients, nodes, metCoeffLag);
fullMatrix<double> *metCoeff;
metCoeff = new fullMatrix<double>(nSampPnts, metCoeffLag.size2());
_bezier->matrixLag2Bez.mult(metCoeffLag, *metCoeff);
@@ -861,33 +816,33 @@ void MetricBasis::_getMetricData(MElement *el, MetricData *&md) const
md = new MetricData(metCoeff, jac, -1, 0, 0);
}
-void MetricBasis::_fillCoeff(const GradientBasis *gradients,
+void MetricBasis::_fillCoeff(int dim, const GradientBasis *gradients,
fullMatrix<double> &nodes, fullMatrix<double> &coeff)
{
const int nSampPnts = gradients->getNumSamplingPoints();
- switch (nodes.size2()) {
+ switch (dim) {
case 0 :
return;
case 1 :
- Msg::Fatal("not implemented");
+ Msg::Fatal("Should not be here, metric for 1d not implemented");
break;
case 2 :
{
- fullMatrix<double> dxydX(nSampPnts,2), dxydY(nSampPnts,2);
+ fullMatrix<double> dxydX(nSampPnts,3), dxydY(nSampPnts,3);
gradients->getGradientsFromNodes(nodes, &dxydX, &dxydY, NULL);
gradients->mapFromIdealElement(&dxydX, &dxydY, NULL);
coeff.resize(nSampPnts, 3);
for (int i = 0; i < nSampPnts; i++) {
- const double &dxdX = dxydX(i,0), &dydX = dxydX(i,1);
- const double &dxdY = dxydY(i,0), &dydY = dxydY(i,1);
- const double dvxdX = dxdX*dxdX + dydX*dydX;
- const double dvxdY = dxdY*dxdY + dydY*dydY;
+ const double &dxdX = dxydX(i,0), &dydX = dxydX(i,1), &dzdX = dxydX(i,2);
+ const double &dxdY = dxydY(i,0), &dydY = dxydY(i,1), &dzdY = dxydY(i,2);
+ const double dvxdX = dxdX*dxdX + dydX*dydX + dzdX*dzdX;
+ const double dvxdY = dxdY*dxdY + dydY*dydY + dzdY*dzdY;
coeff(i, 0) = (dvxdX + dvxdY) / 2;
coeff(i, 1) = dvxdX - coeff(i, 0);
- coeff(i, 2) = (dxdX*dxdY + dydX*dydY);
+ coeff(i, 2) = (dxdX*dxdY + dydX*dydY + dzdX*dzdY);
}
}
break;
@@ -965,72 +920,6 @@ double MetricBasis::_computeMinlagR(const fullVector<double> &jac,
}
}
-int MetricBasis::_paramOnPlane(const fullMatrix<double> &nodes3d,
- fullMatrix<double> &nodes2d)
-{
- // check if in xy-plane
- int i = 0;
- while (i < nodes3d.size1() && nodes3d(i,2) == 0) ++i;
- if (i == nodes3d.size1()) {
- nodes2d.setAsProxy(const_cast<double*>(nodes3d.getDataPtr()),
- nodes3d.size1(), 2);
- return 0;
- }
-
- // check if coplanar and reparameterize
- SVector3 vx(nodes3d(1, 0) - nodes3d(0, 0),
- nodes3d(1, 1) - nodes3d(0, 1),
- nodes3d(1, 2) - nodes3d(0, 2));
- SVector3 vy(nodes3d(2, 0) - nodes3d(0, 0),
- nodes3d(2, 1) - nodes3d(0, 1),
- nodes3d(2, 2) - nodes3d(0, 2));
-
- double lx = vx.norm(), L = vy.norm() + lx;
- if (norm(crossprod(vx, vy)) / L < 1e-10) {
- nodes2d.resize(nodes3d.size1(), 2, true);
- return 0;
- }
-
- nodes2d.resize(nodes3d.size1(), 2, false);
- nodes2d(0, 0) = 0;
- nodes2d(0, 1) = 0;
-
- SVector3 ux = vx;
- ux *= 1/lx;
- nodes2d(1, 0) = lx;
- nodes2d(1, 1) = 0;
-
- double p_vy_x = dot(vy, ux);
- SVector3 uy = ux;
- uy *= -p_vy_x;
- uy += vy;
- double n = uy.norm();
- uy.normalize();
- nodes2d(2, 0) = p_vy_x;
- nodes2d(2, 1) = n;
-
- for (int i = 3; i < nodes3d.size1(); ++i) {
- SVector3 v(nodes3d(i, 0) - nodes3d(0, 0),
- nodes3d(i, 1) - nodes3d(0, 1),
- nodes3d(i, 2) - nodes3d(0, 2));
- double p_v_x = dot(v, ux);
- double p_v_y = dot(v, uy);
- SVector3 px = ux, py = uy;
- px *= -p_v_x;
- py *= -p_v_y;
- v -= px;
- v -= py;
- if (norm(v) / L < 1e-10) {
- //Msg::Info("%g / %g = %g <= 1e-10", norm(v), L, norm(v)/L);
- return 1;
- }
- nodes2d(i, 0) = p_v_x;
- nodes2d(i, 1) = p_v_y;
- }
-
- return 0;
-}
-
void MetricBasis::_minMaxA(
const fullMatrix<double> &coeff, double &min, double &max) const
{
diff --git a/Numeric/MetricBasis.h b/Numeric/MetricBasis.h
index c126aeb3e806de58059ca4675b81c1dac6013950..bab7253e539fe6d8d2975509b8a4ef9be74759bb 100644
--- a/Numeric/MetricBasis.h
+++ b/Numeric/MetricBasis.h
@@ -78,14 +78,11 @@ private:
void _getMetricData(MElement*, MetricData*&) const;
double _subdivideForRmin(MetricData*, double RminLag, double tol) const;
- static void _fillCoeff(const GradientBasis*,
+ static void _fillCoeff(int dim, const GradientBasis*,
fullMatrix<double> &nodes, fullMatrix<double> &coeff);
static double _computeMinlagR(const fullVector<double> &jac,
const fullMatrix<double> &coeff, int num);
- static int _paramOnPlane(const fullMatrix<double> &nodes3d,
- fullMatrix<double> &nodes2d);
-
void _minMaxA(const fullMatrix<double>&, double &min, double &max) const;
void _minK(const fullMatrix<double>&, const fullVector<double>&, double &min) const;
void _maxAstKpos(const fullMatrix<double>&, const fullVector<double>&,
@@ -101,6 +98,7 @@ private:
const double x = .5 * (K - a*a*a + 3*a);
if (x > 1+1e-13 || x < -1-1e-13) {
Msg::Warning("x = %g (|1+%g|)", x, std::abs(x)-1);
+ Msg::Fatal("a");
}
double ans;
diff --git a/Plugin/AnalyseCurvedMesh.cpp b/Plugin/AnalyseCurvedMesh.cpp
index e9671a2d0f82414b06fef4fa68e44df66ee0a896..636575f494caaf4ed863551a269054f4828c54bd 100644
--- a/Plugin/AnalyseCurvedMesh.cpp
+++ b/Plugin/AnalyseCurvedMesh.cpp
@@ -190,9 +190,8 @@ PView *GMSH_AnalyseCurvedMeshPlugin::execute(PView *v)
std::map<int, std::vector<double> > dataPV;
for (unsigned int i = 0; i < _data.size(); ++i) {
MElement *const el = _data[i].element();
- if (el->getDim() == dim) {
+ if (el->getDim() == dim)
dataPV[el->getNum()].push_back(_data[i].minR());
- }
}
std::stringstream name;
name << "min R " << dim << "D";