diff --git a/Numeric/MetricBasis.cpp b/Numeric/MetricBasis.cpp
index 1bc1c72f5c103b842a062c9eef7f247c79325c22..9d293c50908cdbc6eefcb61eeceecb1c5fbe9337 100644
--- a/Numeric/MetricBasis.cpp
+++ b/Numeric/MetricBasis.cpp
@@ -3,6 +3,7 @@
// See the LICENSE.txt file for license information. Please report all
// bugs and problems to the public mailing list <gmsh@geuz.org>.
+#include "AnalyseCurvedMesh.h"
#include "MetricBasis.h"
#include "BasisFactory.h"
#include "pointsGenerators.h"
@@ -53,9 +54,9 @@ MetricBasis::MetricBasis(int tag)
const int metOrder = metricOrder(tag);
if (type == TYPE_HEX || type == TYPE_PRI) {
int order = ElementType::OrderFromTag(tag);
- _jacobian = new JacobianBasis(tag, 3*order);
+ _jacobian = BasisFactory::getJacobianBasis(tag, 3*order);
}
- else if (type == TYPE_TET)
+ else if (type == TYPE_TET || type == TYPE_TRI || type == TYPE_QUA)
_jacobian = BasisFactory::getJacobianBasis(tag);
else
Msg::Fatal("metric not implemented for element tag %d", tag);
@@ -82,19 +83,37 @@ double MetricBasis::minRCorner(MElement *el)
int nSampPnts = jacobian->getNumJacNodes();
if (el->getType() == TYPE_PYR) nSampPnts = 4;
-
int nMapping = gradients->getNumMapNodes();
+
+ // Nodes
fullMatrix<double> nodes(nMapping, 3);
el->getNodesCoord(nodes);
+ fullMatrix<double> nodesMetric;
+ switch (el->getDim()) {
+ case 1:
+ return -1;
- // Metric coefficients
- fullMatrix<double> metCoeffLag;
- _fillCoeff(gradients, nodes, metCoeffLag);
+ 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());
jacobian->getSignedJacobian(nodes, jacLag);
+ // Metric coefficients
+ fullMatrix<double> metCoeffLag;
+ _fillCoeff(gradients, nodesMetric, metCoeffLag);
+
// Compute min_corner(R)
return _computeMinlagR(jacLag, metCoeffLag, nSampPnts);
}
@@ -166,15 +185,27 @@ double MetricBasis::getBoundMinR(MElement *el)
{
int nSampPnts = _gradients->getNumSamplingPoints();
int nMapping = _gradients->getNumMapNodes();
+
+ // Nodes
fullMatrix<double> nodes(nMapping, 3);
el->getNodesCoord(nodes);
+ fullMatrix<double> nodesMetric;
+ switch (el->getDim()) {
+ case 1:
+ return -1;
- // Metric coefficients
- fullMatrix<double> metCoeffLag;
- _fillCoeff(_gradients, nodes, metCoeffLag);
- fullMatrix<double> *metCoeff;
- metCoeff = new fullMatrix<double>(nSampPnts, metCoeffLag.size2());
- _bezier->matrixLag2Bez.mult(metCoeffLag, *metCoeff);
+ 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());
@@ -182,9 +213,16 @@ double MetricBasis::getBoundMinR(MElement *el)
_jacobian->getSignedJacobian(nodes, jacLag);
_jacobian->lag2Bez(jacLag, *jac);
- //
+ // Metric coefficients
+ fullMatrix<double> metCoeffLag;
+ _fillCoeff(_gradients, nodesMetric, metCoeffLag);
+ fullMatrix<double> *metCoeff;
+ metCoeff = new fullMatrix<double>(nSampPnts, metCoeffLag.size2());
+ _bezier->matrixLag2Bez.mult(metCoeffLag, *metCoeff);
+
+ // Compute min(R, _tol)
double RminLag, RminBez;
- _computeRmin(*metCoeff, *jac, RminLag, RminBez, 0);
+ _computeRmin(*metCoeff, *jac, RminLag, RminBez);
if (RminLag-RminBez < MetricBasis::_tol) {
return RminBez;
@@ -288,45 +326,47 @@ void MetricBasis::_fillInequalities(int metricOrder)
}
}
- exp.resize(_jacobian->getBezier()->_exponents.size1(),
- _jacobian->getBezier()->_exponents.size2());
- for (int i = 0; i < _jacobian->getBezier()->_exponents.size1(); ++i) {
- for (int j = 0; j < _jacobian->getBezier()->_exponents.size2(); ++j) {
- exp(i, j) = static_cast<int>(_jacobian->getBezier()->_exponents(i, j) + .5);
+ if (_jacobian) {
+ exp.resize(_jacobian->getBezier()->_exponents.size1(),
+ _jacobian->getBezier()->_exponents.size2());
+ for (int i = 0; i < _jacobian->getBezier()->_exponents.size1(); ++i) {
+ for (int j = 0; j < _jacobian->getBezier()->_exponents.size2(); ++j) {
+ exp(i, j) = static_cast<int>(_jacobian->getBezier()->_exponents(i, j) + .5);
+ }
}
- }
- int njac = _jacobian->getNumJacNodes();
- for (int i = 0; i < njac; i++) {
- for (int j = i; j < njac; j++) {
- int order = metricOrder/2*3;
- double num = 1, den = 1;
- {
- int compl1 = order;
- int compl2 = order;
- int compltot = 2*order;
- for (int k = 0; k < dimSimplex; k++) {
- num *= nChoosek(compl1, exp(i, k))
- * nChoosek(compl2, exp(j, k));
- den *= nChoosek(compltot, exp(i, k) + exp(j, k));
- compl1 -= exp(i, k);
- compl2 -= exp(j, k);
- compltot -= exp(i, k) + exp(j, k);
+ int njac = _jacobian->getNumJacNodes();
+ for (int i = 0; i < njac; i++) {
+ for (int j = i; j < njac; j++) {
+ int order = metricOrder/2*3;
+ double num = 1, den = 1;
+ {
+ int compl1 = order;
+ int compl2 = order;
+ int compltot = 2*order;
+ for (int k = 0; k < dimSimplex; k++) {
+ num *= nChoosek(compl1, exp(i, k))
+ * nChoosek(compl2, exp(j, k));
+ den *= nChoosek(compltot, exp(i, k) + exp(j, k));
+ compl1 -= exp(i, k);
+ compl2 -= exp(j, k);
+ compltot -= exp(i, k) + exp(j, k);
+ }
+ }
+ for (int k = dimSimplex; k < dim; k++) {
+ num *= nChoosek(order, exp(i, k))
+ * nChoosek(order, exp(j, k));
+ den *= nChoosek(2*order, exp(i, k) + exp(j, k));
}
- }
- for (int k = dimSimplex; k < dim; k++) {
- num *= nChoosek(order, exp(i, k))
- * nChoosek(order, exp(j, k));
- den *= nChoosek(2*order, exp(i, k) + exp(j, k));
- }
- if (i != j) num *= 2;
+ if (i != j) num *= 2;
- ++countJ2;
- int hash = 0;
- for (int k = 0; k < dim; k++) {
- hash += (exp(i, k)+exp(j, k)) * pow_int(2*order+1, k);
+ ++countJ2;
+ int hash = 0;
+ for (int k = 0; k < dim; k++) {
+ hash += (exp(i, k)+exp(j, k)) * pow_int(2*order+1, k);
+ }
+ _ineqJ2[hash].push_back(IneqData(num/den, i, j));
}
- _ineqJ2[hash].push_back(IneqData(num/den, i, j));
}
}
@@ -531,8 +571,7 @@ int MetricBasis::metricOrder(int tag)
void MetricBasis::_computeRmin(
const fullMatrix<double> &coeff, const fullVector<double> &jac,
- double &RminLag, double &RminBez,
- int depth, bool debug) const
+ double &RminLag, double &RminBez) const
{
RminLag = _computeMinlagR(jac, coeff, _bezier->getNumLagCoeff());
if (RminLag == 0) {
@@ -540,6 +579,13 @@ void MetricBasis::_computeRmin(
return;
}
+ if (coeff.size2() == 3) { // 2d element
+ double mina, dummy;
+ _minMaxA(coeff, mina, dummy);
+ RminBez = std::sqrt(_R2Dsafe(mina));
+ return;
+ }
+
double minK;
_minK(coeff, jac, minK);
if (minK < 1e-10) {
@@ -595,9 +641,9 @@ void MetricBasis::_computeRmin(
double am1 = std::pow(amApprox*amApprox*amApprox+da+.05, 1/3.);
//double am0S = am0, am1S = am1;
double am = (am0 + am1)/2;
- double R0 = _Rsafe(am0, minK);
- double R1 = _Rsafe(am1, minK);
- double Rnew = _Rsafe(am, minK);
+ double R0 = _R3Dsafe(am0, minK);
+ double R1 = _R3Dsafe(am1, minK);
+ double Rnew = _R3Dsafe(am, minK);
int cnt = 0;
while (std::abs(R0-Rnew) > _tol*.01 || std::abs(R1-Rnew) > _tol*.01) {
@@ -611,17 +657,17 @@ void MetricBasis::_computeRmin(
R1 = Rnew;
}
am = (am0 + am1)/2;
- Rnew = _Rsafe(am, minK);
+ Rnew = _R3Dsafe(am, minK);
}
if (am < maxa) {
- RminBez = _Rsafe(am, minK);
+ RminBez = _R3Dsafe(am, minK);
RminBez = std::sqrt(RminBez);
return;
}
}
- RminBez = _Rsafe(maxa, minK);
+ RminBez = _R3Dsafe(maxa, minK);
RminBez = std::sqrt(RminBez);
return;
}
@@ -670,7 +716,7 @@ void MetricBasis::_computeRmax(
RmaxLag = std::max(RmaxLag, std::sqrt((a - std::sqrt(3)) / (a + std::sqrt(3))));
}
else {
- const double tmpR = _Rsafe(a, jac(i)/p/p*jac(i)/p);
+ const double tmpR = _R3Dsafe(a, jac(i)/p/p*jac(i)/p);
RmaxLag = std::max(RmaxLag, std::sqrt(tmpR));
}
}
@@ -707,7 +753,7 @@ double MetricBasis::_subdivideForRmin(MetricData *md, double RminLag, double tol
jac = new fullVector<double>;
jac->setAsProxy(*subjac, i * numJacPnts, numJacPnts);
double minLag, minBez;
- _computeRmin(*coeff, *jac, minLag, minBez, depth+1);
+ _computeRmin(*coeff, *jac, minLag, minBez);
RminLag = std::min(RminLag, minLag);
int newNum = num + (i+1) * pow_int(10, depth);
MetricData *metData = new MetricData(coeff, jac, minBez, depth+1, newNum);
@@ -771,15 +817,31 @@ void MetricBasis::_getMetricData(MElement *el, MetricData *&md) const
{
int nSampPnts = _gradients->getNumSamplingPoints();
int nMapping = _gradients->getNumMapNodes();
+
+ // Nodes
fullMatrix<double> nodes(nMapping, 3);
el->getNodesCoord(nodes);
+ fullMatrix<double> nodesMetric;
+ switch (el->getDim()) {
+ case 1:
+ md = NULL;
+ return;
- // Metric coefficients
- fullMatrix<double> metCoeffLag;
- _fillCoeff(_gradients, nodes, metCoeffLag);
- fullMatrix<double> *metCoeff;
- metCoeff = new fullMatrix<double>(nSampPnts, metCoeffLag.size2());
- _bezier->matrixLag2Bez.mult(metCoeffLag, *metCoeff);
+ 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());
@@ -787,6 +849,13 @@ void MetricBasis::_getMetricData(MElement *el, MetricData *&md) const
_jacobian->getSignedJacobian(nodes, jacLag);
_jacobian->lag2Bez(jacLag, *jac);
+ // Metric coefficients
+ fullMatrix<double> metCoeffLag;
+ _fillCoeff(_gradients, nodesMetric, metCoeffLag);
+ fullMatrix<double> *metCoeff;
+ metCoeff = new fullMatrix<double>(nSampPnts, metCoeffLag.size2());
+ _bezier->matrixLag2Bez.mult(metCoeffLag, *metCoeff);
+
md = new MetricData(metCoeff, jac, -1, 0, 0);
}
@@ -799,10 +868,27 @@ void MetricBasis::_fillCoeff(const GradientBasis *gradients,
case 0 :
return;
case 1 :
- case 2 :
Msg::Fatal("not implemented");
break;
+ case 2 :
+ {
+ fullMatrix<double> dxydX(nSampPnts,2), dxydY(nSampPnts,2);
+ gradients->getGradientsFromNodes(nodes, &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;
+ coeff(i, 0) = (dvxdX + dvxdY) / 2;
+ coeff(i, 1) = dvxdX - coeff(i, 0);
+ coeff(i, 2) = (dxdX*dxdY + dydX*dydY);
+ }
+ }
+ break;
+
case 3 :
{
fullMatrix<double> dxyzdX(nSampPnts,3), dxyzdY(nSampPnts,3), dxyzdZ(nSampPnts,3);
@@ -834,28 +920,112 @@ double MetricBasis::_computeMinlagR(const fullVector<double> &jac,
const fullMatrix<double> &coeff, int num)
{
double Rmin = 1.;
- for (int i = 0; i < num; ++i) {
- if (jac(i) <= 0.) {
- return 0;
- }
- else {
- double q = coeff(i, 0);
+
+ switch (coeff.size2()) {
+ case 7:
+ for (int i = 0; i < num; ++i) {
+ if (jac(i) <= 0.) return 0;
+
+ const double q = coeff(i, 0);
double p = 0;
for (int k = 1; k < 7; ++k) {
p += pow_int(coeff(i, k), 2);
}
p = std::sqrt(p/6);
const double a = q/p;
- if (a > 1e4) {
+ if (a > 1e4) { // TODO: from _tol ?
Rmin = std::min(Rmin, std::sqrt((a - std::sqrt(3)) / (a + std::sqrt(3))));
}
else {
- const double tmpR = _Rsafe(a, jac(i)/p/p*jac(i)/p);
+ const double tmpR = _R3Dsafe(a, jac(i)/p/p*jac(i)/p);
Rmin = std::min(Rmin, std::sqrt(tmpR));
}
}
+ return Rmin;
+
+ case 3:
+ for (int i = 0; i < num; ++i) {
+ if (jac(i) <= 0.) return 0;
+
+ const double q = coeff(i, 0);
+ const double p = pow_int(coeff(i, 1), 2) + pow_int(coeff(i, 2), 2);
+ const double a = q/std::sqrt(p);
+ const double tmpR = _R2Dsafe(a);
+ Rmin = std::min(Rmin, std::sqrt(tmpR));
+ }
+ return Rmin;
+
+ default:
+ Msg::Error("coeff have not right number of column");
+ return -1;
+ }
+}
+
+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);
+ Msg::Info("Was on xy-plane");
+ 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;
}
- return Rmin;
+
+ 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(
@@ -871,7 +1041,7 @@ void MetricBasis::_minMaxA(
const int i = it->second[k].i;
const int j = it->second[k].j;
double tmp = 0;
- for (int l = 1; l < 7; ++l) {
+ for (int l = 1; l < coeff.size2(); ++l) {
tmp += coeff(i, l) * coeff(j, l);
}
den += it->second[k].val * tmp;
@@ -882,8 +1052,10 @@ void MetricBasis::_minMaxA(
max = std::max(val, max);
++it;
}
- min *= 6;
- max *= 6;
+ if (coeff.size2() == 7) {
+ min *= 6;
+ max *= 6;
+ }
min = min > 1 ? std::sqrt(min) : 1;
max = std::sqrt(max);
diff --git a/Numeric/MetricBasis.h b/Numeric/MetricBasis.h
index 8bbe4c4aa049a0b6e29e47476649282ca7b78c67..c126aeb3e806de58059ca4675b81c1dac6013950 100644
--- a/Numeric/MetricBasis.h
+++ b/Numeric/MetricBasis.h
@@ -70,7 +70,7 @@ private:
void _lightenInequalities(int&, int&, int&); //TODO change
void _computeRmin(const fullMatrix<double>&, const fullVector<double>&,
- double &RminLag, double &RminBez, int depth, bool debug = false) const;
+ double &RminLag, double &RminBez) const;
void _computeRmax(const fullMatrix<double>&, const fullVector<double>&,
double &RmaxLag) const;
void _computeTermBeta(double &a, double &K, double &dRda,
@@ -83,6 +83,9 @@ private:
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>&,
@@ -94,7 +97,7 @@ private:
void _maxKstAsharp(const fullMatrix<double>&, const fullVector<double>&,
double mina, double beta, double &maxK) const;
- static double _Rsafe(double a, double K) {
+ static double _R3Dsafe(double a, double K) {
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);
@@ -115,6 +118,10 @@ private:
}
return ans;
}
+ static double _R2Dsafe(double a) {
+ if (a < 1) Msg::Warning("R2d = %g", (a - 1) / (a + 1));
+ return (a - 1) / (a + 1);
+ }
private:
class gterIneq {