diff --git a/Geo/MElement.cpp b/Geo/MElement.cpp
index aca031d8ee093733809b0cefe3861ede0df9cf61..d6322bdacb311cb8233437626ed0847e79b332bd 100644
--- a/Geo/MElement.cpp
+++ b/Geo/MElement.cpp
@@ -235,10 +235,10 @@ double MElement::maxDistToStraight() const
return maxdx;
}
-double MElement::minAnisotropyMeasure(bool knownValid, bool reversedOK)
+double MElement::minIsotropyMeasure(bool knownValid, bool reversedOK)
{
#if defined(HAVE_MESH)
- return jacobianBasedQuality::minAnisotropyMeasure(this, knownValid, reversedOK);
+ return jacobianBasedQuality::minIsotropyMeasure(this, knownValid, reversedOK);
#else
return 0.;
#endif
diff --git a/Geo/MElement.h b/Geo/MElement.h
index 32bd4a1851509814d5dd873029af9395a2a11921..8c3b84b5e045d323d17b65017124d035a88bccf1 100644
--- a/Geo/MElement.h
+++ b/Geo/MElement.h
@@ -218,7 +218,7 @@ class MElement
signedInvCondNumRange(sICNMin, sICNMax);
return sICNMin;
}
- double minAnisotropyMeasure(bool knownValid = false, bool reversedOk = false);
+ double minIsotropyMeasure(bool knownValid = false, bool reversedOk = false);
double minScaledJacobian(bool knownValid = false, bool reversedOk = false);
double specialQuality();
double specialQuality2();
diff --git a/Mesh/qualityMeasuresJacobian.cpp b/Mesh/qualityMeasuresJacobian.cpp
index bae81043ca4f4ae67e637841d6f4d8cc9cfd23f1..81c82855b25d85d7248832da8b75f9528dd32ce8 100644
--- a/Mesh/qualityMeasuresJacobian.cpp
+++ b/Mesh/qualityMeasuresJacobian.cpp
@@ -27,12 +27,6 @@ double TM##n = Cpu();
namespace jacobianBasedQuality {
-bool _CoeffDataAnisotropy::hasError = false;
-MElement *_CoeffDataAnisotropy::currentElement = NULL;
-_CoeffDataAnisotropy *_CoeffDataAnisotropy::current = NULL;
-std::fstream _CoeffDataAnisotropy::file;
-std::vector<double> _CoeffDataAnisotropy::mytimes;
-std::vector<int> _CoeffDataAnisotropy::mynumbers;
double _CoeffDataScaledJac::cTri = 2/std::sqrt(3);
double _CoeffDataScaledJac::cTet = std::sqrt(2);
double _CoeffDataScaledJac::cPyr = std::sqrt(2);
@@ -40,7 +34,6 @@ double _CoeffDataScaledJac::cPyr = std::sqrt(2);
void minMaxJacobianDeterminant(MElement *el, double &min, double &max,
const fullMatrix<double> *normals)
{
- _CoeffDataAnisotropy::currentElement = el;
const JacobianBasis *jfs = el->getJacobianFuncSpace();
if (!jfs) {
Msg::Error("Jacobian function space not implemented for type of element %d", el->getTypeForMSH());
@@ -171,110 +164,6 @@ double minScaledJacobian(MElement *el, bool knownValid, bool reversedOk)
return min;
}
-double minAnisotropyMeasure(MElement *el,
- bool knownValid,
- bool reversedOk,
- int n)
-{
- if (!knownValid) {
- double jmin, jmax;
- minMaxJacobianDeterminant(el, jmin, jmax);
- if (jmin <= 0 && jmax >= 0) return 0;
- if (!reversedOk && jmax < 0) return 0;
- }
- // Computation of the minimum of the anisotropy measure should never
- // be performed to invalid elements (for which the measure is 0).
-
- if (_CoeffDataAnisotropy::noMoreComputationPlease()) {
- Msg::Info("Sorry, no more computation");
- return -9; //fordebug
- }
-
-// double minSampled = minSampledAnisotropyMeasure(el, n); //fordebug
-
- fullMatrix<double> nodesXYZ(el->getNumVertices(), 3);
- el->getNodesCoord(nodesXYZ);
-
- const GradientBasis *gradBasis;
- const JacobianBasis *jacBasis = NULL;
-
- const int type = el->getType();
- const bool serendipFalse = false;
- const int metricOrder = _CoeffDataAnisotropy::metricOrder(el);
- const int jacDetOrder = 3*metricOrder/2;
- if (metricOrder == -1) {
- Msg::Error("Anisotropy measure not implemented for type of element %d", el->getType());
- return -1;
- }
-
- FuncSpaceData metricSpace, jacDetSpace;
-
- switch(type) {
- case TYPE_TET:
- case TYPE_HEX:
- case TYPE_PRI:
- jacDetSpace = FuncSpaceData(el, jacDetOrder, &serendipFalse);
- jacBasis = BasisFactory::getJacobianBasis(jacDetSpace);
- //no break
- case TYPE_TRI:
- case TYPE_QUA:
- metricSpace = FuncSpaceData(el, metricOrder, &serendipFalse);
- break;
- case TYPE_PYR:
- metricSpace = FuncSpaceData(el, false, metricOrder+2, metricOrder, &serendipFalse);
- jacDetSpace = FuncSpaceData(el, false, jacDetOrder+3, jacDetOrder, &serendipFalse);
- jacBasis = BasisFactory::getJacobianBasis(jacDetSpace);
- break;
- }
- gradBasis = BasisFactory::getGradientBasis(metricSpace);
-
- // Metric Coefficients
- fullMatrix<double> coeffMetricBez;
- {
- fullMatrix<double> coeffMetricLag;
- _CoeffDataAnisotropy::fillMetricCoeff(gradBasis, nodesXYZ, coeffMetricLag,
- el->getDim(), true);
-
- coeffMetricBez.resize(coeffMetricLag.size1(), coeffMetricLag.size2());
- gradBasis->lag2Bez(coeffMetricLag, coeffMetricBez);
- }
-
- // Jacobian determinant coefficients
- fullVector<double> coeffJacDetBez;
- if (jacBasis) {
- fullVector<double> coeffDetLag(jacBasis->getNumJacNodes());
- jacBasis->getSignedIdealJacobian(nodesXYZ, coeffDetLag);
-
- coeffJacDetBez.resize(jacBasis->getNumJacNodes());
- jacBasis->lag2Bez(coeffDetLag, coeffJacDetBez);
- }
-
- std::vector<_CoeffData*> domains;
- domains.push_back(
- new _CoeffDataAnisotropy(coeffJacDetBez, coeffMetricBez,
- jacBasis ? jacBasis->getBezier() : NULL,
- gradBasis->getBezier(), 0, el)
- );
-
- _subdivideDomains(domains);
- if (domains.size()/7 > 500) {//fordebug
- Msg::Info("A too much subdivision: %d (el %d, type %d)", domains.size()/7, el->getNum(), el->getType());
- }
-
- double min = domains[0]->minB();
- delete domains[0];
- for (unsigned int i = 1; i < domains.size(); ++i) {
- min = std::min(min, domains[i]->minB());
- delete domains[i];
- }
-
-// if (min - minSampled > 1e-13) {//fordebug
-// Msg::Error("no, works no el %d (%g vs %g, diff %g)", el->getNum(), min, minSampled, min-minSampled);
-// }
-
- return min;
-}
-
double minIsotropyMeasure(MElement *el,
bool knownValid,
bool reversedOk)
@@ -1063,809 +952,6 @@ double _CoeffDataIsotropy::_computeLowerBound() const
// coeffDenominator, true), 1/3.);
}
-// Anisotropy measure
-_CoeffDataAnisotropy::_CoeffDataAnisotropy(fullVector<double> &det,
- fullMatrix<double> &metric,
- const bezierBasis *bfsDet,
- const bezierBasis *bfsMet,
- int depth,
- MElement *el,
- int num)
-: _CoeffData(depth), _coeffsJacDet(det.getDataPtr(), det.size()),
- _coeffsMetric(metric.getDataPtr(), metric.size1(), metric.size2()),
- _bfsDet(bfsDet), _bfsMet(bfsMet), _elForDebug(el), _numForDebug(num)
-{
- _CoeffDataAnisotropy::current = this;
- if (!det.getOwnData() || !metric.getOwnData()) {
- Msg::Fatal("Cannot create an instance of _CoeffDataScaledJac from a "
- "fullVector or a fullMatrix that does not own its data.");
- }
- // _coeffsJacDet and _coeffsMetric reuse data, this avoid to allocate new
- // arrays and to copy data that are not used outside of this object.
- // It remains to swap ownerships:
- det.setOwnData(false);
- metric.setOwnData(false);
- const_cast<fullVector<double>&>(_coeffsJacDet).setOwnData(true);
- const_cast<fullMatrix<double>&>(_coeffsMetric).setOwnData(true);
-
- _computeAtCorner(_minL, _maxL);
-
- _minB = 0;
- if (boundsOk(_minL, _maxL)) return;
- else _minB = _computeLowerBound();
-
- // computation of _maxB not implemented for now
- //statsForMatlab(_elForDebug, _numForDebug);//fordebug
-}
-
-void _CoeffDataAnisotropy::_computeAtCorner(double &min, double &max) const
-{
- min = std::numeric_limits<double>::infinity();
- max = -min;
-
- for (int i = 0; i < _bfsMet->getNumLagCoeff(); i++) {
- const double &q = _coeffsMetric(i, 0);
- double p = 0;
- for (int k = 1; k < _coeffsMetric.size2(); ++k) {
- p += pow_int(_coeffsMetric(i, k), 2);
- }
- p = std::sqrt(p);
- double qualSqrd;
- if (_coeffsMetric.size2() == 3)
- qualSqrd = _R2Dsafe(q, p);
- else
- qualSqrd = _R3Dsafe(q, p, _coeffsJacDet(i));
- min = std::min(min, qualSqrd);
- max = std::max(max, qualSqrd);
- }
- min = std::sqrt(min);
- max = std::sqrt(max);
-}
-
-double _CoeffDataAnisotropy::_computeLowerBound() const
-{
- // 2D element
- if (_coeffsMetric.size2() == 3) {
- double mina = _computeLowerBoundA();
- return std::sqrt(_R2Dsafe(mina));
- }
-
- // 3D element
- double mina = _computeLowerBoundA();
- double minK = _computeLowerBoundK();
-
- _moveInsideDomain(mina, minK, true);
-
- double p_dRda, p_dRdK;
- _computePseudoDerivatives(mina, minK, p_dRda, p_dRdK);
-
- if (p_dRda < 0) {
- // cases 1 & 2 => compute a lower bounding curve K = beta * a^3 + c * a
- // with c such that the curve that passes through (minK, mina) has the
- // slope equal to -dRda/dRdK
- double slope = -p_dRda/p_dRdK;
- double gamma = .5*(3*minK/mina - slope);
- double beta;
- _computeBoundingCurve(beta, gamma, true);
- /*{//fordebug
- double beta = (minK/mina-c)/mina/mina;
- if (std::abs(p_dRda + p_dRdK * (3*beta*mina*mina+c)) > 1e-12) {
- Msg::Error("Derivative along curve not zero %g", p_dRda + p_dRdK * (3*beta*mina*mina+c));
- }
- }*/
-
- double a_int = mina, K_int = minK;
- int where = _intersectionCurveLeftCorner(beta, gamma, a_int, K_int);
- if (where == 1) {
- // the minimum is on the curve, long to compute, we return this for now:
- return std::sqrt(_R3Dsafe(mina, minK));
- }
- //Msg::Info("int (%g, %g) (beta %g, gamma %g)", a_int, K_int, beta, gamma);//fordebug
-
- // Let a0 be the point at which R(a, minK) takes its minimum. If there is
- // an intersection and if a_int < a0, the minimum is at (a_int, minK),
- // otherwise it is at (a0, minK).
- bool minIsAta0;
- if (where == -1 || _moveInsideDomain(a_int, K_int, false)) {
- minIsAta0 = true;
- }
- else {
- double p_dRda, p_dRdK;
- _computePseudoDerivatives(a_int, K_int, p_dRda, p_dRdK);
- /*if (p_dRda + p_dRdK * (3*beta*a_int*a_int+c) < -1e-5) {//fordebug
- Msg::Error("Derivative along curve not positive or zero %g", p_dRda + p_dRdK * (3*beta*a_int*a_int+c));
- Msg::Info("(%g, %g) vs (%g, %g), diff (%g, %g)", mina, minK, a_int, K_int, a_int-mina, K_int-minK);
- double beta2 = (minK/mina-c)/mina/mina;
- Msg::Info("beta %g - %g = %g", beta, beta2, beta-beta2);
- }*/
- minIsAta0 = p_dRda > 0;
- }
-
- if (minIsAta0) {
- double a0 = _computeAbscissaMinR(mina, minK);
- return std::sqrt(_R3Dsafe(a0, minK));
- }
- else {
- return std::sqrt(_R3Dsafe(a_int, minK));
- }
- }
- else if (p_dRdK < 0) {
- // cases 4 & 5 => compute an upper bounding curve K = beta * a^3 + c * a
- // with c such that the curve that passes through (minK, mina) has the
- // slope equal to -dRda/dRdK
- double slope = -p_dRda/p_dRdK;
- double gamma = .5*(3*minK/mina - slope);
- double beta;
- _computeBoundingCurve(beta, gamma, false);
-
- double a_int = mina, K_int = minK;
- int where = _intersectionCurveLeftCorner(beta, gamma, a_int, K_int);
- if (where == 0) {
- // the minimum is on the curve, long to compute, we return this for now:
- return std::sqrt(_R3Dsafe(mina, minK));
- }
-
- // Let K0 be the point at which R(mina, K) takes its minimum. If there is
- // an intersection and if K_int < K0, the minimum is at (mina, K_int),
- // otherwise it is at (mina, K0).
- double K0 = 2*std::cos(3*std::acos(-1/mina)-M_PI) + (mina*mina-3) * mina;
- if (where == -1)
- return std::sqrt(_R3Dsafe(mina, K0));
- else
- return std::sqrt(_R3Dsafe(mina, std::min(K0, K_int)));
- }
- else {
- return std::sqrt(_R3Dsafe(mina, minK));
- }
-}
-
-double _CoeffDataAnisotropy::_computeLowerBoundA() const
-{
- fullVector<double> coeffNumerator;
- {
- fullVector<double> coeffq;
- coeffq.setAsProxy(_coeffsMetric, 0);
- _bfsMet->getRaiser()->computeCoeff(coeffq, coeffq, coeffNumerator);
- }
-
- fullVector<double> coeffDenominator;
- {
- fullMatrix<double> prox, terms;
- prox.setAsProxy(_coeffsMetric, 1, _coeffsMetric.size2()-1);
- _bfsMet->getRaiser()->computeCoeff(prox, prox, terms);
- coeffDenominator.resize(terms.size1(), true);
- for (int i = 0; i < terms.size1(); ++i) {
- for (int j = 0; j < terms.size2(); ++j) {
- coeffDenominator(i) += terms(i, j);
- }
- }
- }
-
- double bound = _computeBoundRational(coeffNumerator, coeffDenominator, true);
- return bound > 1 ? std::sqrt(bound) : 1;
-}
-
-double _CoeffDataAnisotropy::_computeLowerBoundK() const
-{
- fullVector<double> coeffNumerator;
- _bfsDet->getRaiser()->computeCoeff(_coeffsJacDet,
- _coeffsJacDet,
- coeffNumerator);
-
- fullVector<double> coeffDenominator;
- {
- fullVector<double> P(_coeffsMetric.size1());
- // element of P are automatically set to 0
- for (int i = 0; i < _coeffsMetric.size1(); ++i) {
- for (int l = 1; l < 7; ++l) {
- P(i) += _coeffsMetric(i, l) * _coeffsMetric(i, l);
- }
- P(i) = std::sqrt(P(i));
- }
- _bfsMet->getRaiser()->computeCoeff(P, P, P, coeffDenominator);
- }
-
- return _computeBoundRational(coeffNumerator, coeffDenominator, true);
-}
-
-void _CoeffDataAnisotropy::_computeBoundingCurve(double &beta, double gamma,
- bool compLowBound) const
-{
- // compute a lower/upper bounding curve K = beta * a^3 + gamma * a
- // with gamma fixed.
-
- fullVector<double> J2;
- _bfsDet->getRaiser()->computeCoeff(_coeffsJacDet, _coeffsJacDet, J2);
-
- fullVector<double> q(_coeffsMetric.size1());
- for (int i = 0; i < q.size(); ++i) {
- q(i) = _coeffsMetric(i, 0);
- }
-
- fullVector<double> qp2;
- {
- fullMatrix<double> terms_p, terms_qp2;
- terms_p.setAsProxy(_coeffsMetric, 1, _coeffsMetric.size2()-1);
- _bfsMet->getRaiser()->computeCoeff(q, terms_p, terms_p, terms_qp2);
- qp2.resize(terms_qp2.size1(), true);
- for (int i = 0; i < terms_qp2.size1(); ++i) {
- for (int j = 0; j < terms_qp2.size2(); ++j) {
- qp2(i) += terms_qp2(i, j);
- }
- }
- }
-
- fullVector<double> &coeffNumerator = J2;
- coeffNumerator.axpy(qp2, -gamma);
-
- fullVector<double> coeffDenominator;
- _bfsMet->getRaiser()->computeCoeff(q, q, q, coeffDenominator);
-
- beta = _computeBoundRational(coeffNumerator, coeffDenominator, compLowBound);
-}
-
-void _CoeffDataAnisotropy::fillMetricCoeff(const GradientBasis *gradients,
- const fullMatrix<double> &nodes,
- fullMatrix<double> &coeff,
- int dim, bool ideal)
-{
- const int nSampPnts = gradients->getNumSamplingPoints();
-
- switch (dim) {
- case 0 :
- return;
-
- case 1 :
- Msg::Fatal("Should not be here, metric for 1d not implemented");
- break;
-
- case 2 :
- {
- fullMatrix<double> dxydX(nSampPnts,3);
- fullMatrix<double> dxydY(nSampPnts,3);
- if (ideal)
- gradients->getIdealGradientsFromNodes(nodes, &dxydX, &dxydY, NULL);
- else
- 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);
- double dzdX = 0, dzdY = 0;
- if (nodes.size2() > 2) {
- dzdX = dxydX(i,2);
- 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 + dzdX*dzdY);
- }
- }
- break;
-
- case 3 :
- {
- fullMatrix<double> dxyzdX(nSampPnts,3);
- fullMatrix<double> dxyzdY(nSampPnts,3);
- fullMatrix<double> dxyzdZ(nSampPnts,3);
- if (ideal)
- gradients->getIdealGradientsFromNodes(nodes, &dxyzdX, &dxyzdY, &dxyzdZ);
- else
- gradients->getGradientsFromNodes(nodes, &dxyzdX, &dxyzdY, &dxyzdZ);
-
- coeff.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);
- 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;
- coeff(i, 0) = (dvxdX + dvxdY + dvxdZ) / 3;
- static double fact1 = 1./std::sqrt(6.);
- static double fact2 = 1./std::sqrt(3.);
- coeff(i, 1) = fact1 * (dvxdX - coeff(i, 0));
- coeff(i, 2) = fact1 * (dvxdY - coeff(i, 0));
- coeff(i, 3) = fact1 * (dvxdZ - coeff(i, 0));
- coeff(i, 4) = fact2 * (dxdX*dxdY + dydX*dydY + dzdX*dzdY);
- coeff(i, 5) = fact2 * (dxdZ*dxdY + dydZ*dydY + dzdZ*dzdY);
- coeff(i, 6) = fact2 * (dxdX*dxdZ + dydX*dydZ + dzdX*dzdZ);
- }
- }
- break;
- }
-}
-
-bool _CoeffDataAnisotropy::_moveInsideDomain(double &a,
- double &K,
- bool bottomleftCorner)
-{
- // - bottomleftCorner: if false, then it is a top left or a bottom right corner
- // (it is the same treatment)
- // - Return true if the point has been moved, false otherwise
-
- // NB: When moving 'a', we use N-R to compute the root of
- // f(a) = .5 * (K - a^3 + 3*a) - 'w' where 'w' is -1 or 1.
-
- if (bottomleftCorner) {
- K = std::max(K, .0);
- a = std::max(a, 1.);
- }
-
- double w = _w(a, K);
- if (std::abs(w) <= 1) return false;
-
- const double tolerance = 1e-5;
-
- if (w > 1) {
- if (bottomleftCorner) {
- // make sure to compute the right root (a>1) and
- // avoid slopes that are too small:
- a = std::max(a, 1.1);
-
- // N-R approaches the root from the right, we can impose w > 1
- while (w > 1 || w < 1-tolerance) {
- double df = w - 1;
- double slope = 1.5*(1-a*a);
- a -= df/slope;
- w = _w(a, K);
- }
- }
- else {
- K -= 2*w - 2;
- }
- }
- else {
- if (bottomleftCorner) {
- K -= 2*w + 2;
- }
- else {
- // make sure to compute the right root (a>2)
- a = std::max(a, 2.);
-
- while (std::abs(w+1) > tolerance) {
- double df = w + 1;
- double slope = 1.5*(1-a*a);
- a -= df/slope;
- w = _w(a, K);
- }
- }
- }
- return true;
-}
-
-bool _CoeffDataAnisotropy::_computePseudoDerivatives(double a, double K,
- double &dRda, double &dRdK)
-{
- // Return derivative without the positive (but non-constant) coefficient
- // Useful when only the sign of dRda and dRdK or the ratio dRda/dRdK
- // are needed.
-
- double w = _wSafe(a, K);
- if (std::abs(w) > 1) {
- Msg::Error("Cannot compute pseudo derivatives of R "
- "outside the domain (w=%g)", w);
- return false;
- }
-
- const double phip = (std::acos(w) + M_PI) / 3;
- const double C = 1 + a * std::cos(phip);
-
- dRdK = C / 3;
- dRda = 2 * std::sqrt(1-w*w) * std::sin(phip) + (1-a*a) * C;
- return true;
-}
-
-bool _CoeffDataAnisotropy::_computeDerivatives(double a, double K,
- double &dRda, double &dRdK,
- double &dRdaa, double &dRdaK,
- double &dRdKK)
-{
- const double w = _wSafe(a, K);
- if (std::abs(w) > 1) {
- Msg::Error("Cannot compute derivatives of R outside the domain (w=%g)", w);
- return false;
- }
-
- const double phi = (std::acos(w)) / 3;
- const double phip = phi + M_PI / 3;
- const double S = 1./std::sqrt(1-w*w);
- const double A = (1-a*a)*S;
- const double sin = std::sin(phi);
- const double sinp = std::sin(phip);
- const double cosp = std::cos(phip);
- const double C = 1 + a*cosp;
- const double D = 1./(a+2*std::cos(phi));
- static const double sq3 = std::sqrt(3);
- const double coeff = sq3*D*D/18;
-
- dRdK = coeff*6 * (C*S);
- dRda = coeff*18 * (2*sinp + A*C);
- dRdKK = coeff*S*S * (a*sinp - 4*C*sin*D + 3*w*C*S);
- dRdaK = coeff*S*3 * (a*A*sinp + 3*w*A*C*S + 2*cosp - 4*C*(1+A*sin)*D);
- dRdaa = coeff*9 * (3*w*A*A*C*S - 4*a*C*S + a*A*A*sinp - 4*(1+A*sin)*D*(2*sinp + A*C));
- return true;
-}
-
-int _CoeffDataAnisotropy::_intersectionCurveLeftCorner(double beta, double gamma,
- double &a, double &K)
-{
- // curve K = beta * a^3 + gamma * a
- // on input, a and K are the bottom left corner
- // on output, a and K are the intersection
- // return 0 if the intersection is on the horizontal
- // 1 if the intersection is on the vertical
- // -1 if there is no intersection
-
-// if (beta < 0) {
-// Msg::Warning("Negative curvature, there are 2 or 0 intersections... Is it "
-// "normal? Anyway, I will compute the left one if it exists.");
-//// Msg::Info("%g %g | %g %g", beta, gamma, a, K);
-// }
-
- const double minK = K;
- K = (beta*a*a + gamma)*a;
- if (K >= minK) return 1;
-
- // Find the root by Newton-Raphson.
- K = minK;
- if (beta < 0) {
- // When beta < 0, check if the intersection exists:
- double aMaximum = std::sqrt(-gamma/3/beta);
- double KMaximum = (beta * aMaximum*aMaximum + gamma) * aMaximum;
- if (gamma < 0 || aMaximum < a || KMaximum < K) {
-// Msg::Warning("Sorry but there is no intersection");
- return -1;
- }
- }
- else if (beta > 0) {
- // When beta > 0 and aMin > 'a', we must move 'a' to the right of aMin,
- // otherwise we would compute a wrong intersection:
- double aMinimum = std::sqrt(-gamma/3/beta);
- if (aMinimum > a) a += 2*(aMinimum - a);
- }
-
- // Which precision? We know that w is in [-1, 1] with w = .5 * (K + (3-a*a)*a)
- // We thus seek a good precision on 'K' and '3*a - a^3'
- double a3 = (3-a*a)*a;
- double da3 = std::numeric_limits<double>::infinity();
- double dK = std::numeric_limits<double>::infinity();
- while (std::abs(da3) > 1e-5 || std::abs(dK) > 1e-4) {
- const double slope = 3*beta*a*a + gamma;
- dK = (beta*a*a + gamma)*a - minK;
- a -= dK / slope;
- da3 = a3;
- a3 = (3-a*a)*a;
- da3 -= a3;
- }
- return 0;
-}
-
-double _CoeffDataAnisotropy::_computeAbscissaMinR(double aStart, double K)
-{
- // If K == 0, then R == 0. Note, there are less numerical problems when
- // computing R(2, 0) than R(1, 0), so return 2:
- if (K == 0) return 2;
-
- double a = aStart;
- double a3 = (3-a*a)*a;
- double da3 = std::numeric_limits<double>::infinity();
- double dRda = da3, dRdK, dRdaa, dRdaK, dRdKK;
- while (std::abs(da3) > 1e-5 || std::abs(dRda) > 1e-5) {
- _computeDerivatives(a, K, dRda, dRdK, dRdaa, dRdaK, dRdKK);
- double da = - dRda / dRdaa;
- da3 = a3;
- double potentiala = a+da;
- a3 = (3-potentiala*potentiala)*potentiala;
- da3 -= a3;
- // Make sure we are not going outside the domain:
- while (std::abs(_w(a+da, K)) > 1) da /= 2;
- a += da;
- }
- return a;
-}
-
-bool _CoeffDataAnisotropy::boundsOk(double minL, double maxL) const
-{
- static double tol = 1e-3;
- return minL < tol*1e-3 || (_minB > minL-tol && _minB > minL*(1-100*tol));
-}
-
-void _CoeffDataAnisotropy::getSubCoeff(std::vector<_CoeffData*> &v) const
-{
- v.clear();
- v.reserve(_bfsMet->getNumDivision());
- fullMatrix<double> subCoeffM;
- fullVector<double> subCoeffD;
- _bfsMet->subdivideBezCoeff(_coeffsMetric, subCoeffM);
- // in 2D, _bfsDet is NULL
- if (_bfsDet) _bfsDet->subdivideBezCoeff(_coeffsJacDet, subCoeffD);
-
- // in 2D, _coeffsJacDet is empty (szD == 0)
- int szD = _coeffsJacDet.size();
- int szM1 = _coeffsMetric.size1();
- int szM2 = _coeffsMetric.size2();
- for (int i = 0; i < _bfsMet->getNumDivision(); i++) {
- fullVector<double> coeffD(szD);
- fullMatrix<double> coeffM(szM1, szM2);
- coeffD.copy(subCoeffD, i * szD, szD, 0);
- coeffM.copy(subCoeffM, i * szM1, szM1, 0, szM2, 0, 0);
- int newNum = _numForDebug + (i+1) * pow_int(10, _depth);
- _CoeffDataAnisotropy *newData
- = new _CoeffDataAnisotropy(coeffD, coeffM, _bfsDet, _bfsMet, _depth+1, _elForDebug, newNum);
- v.push_back(newData);
- }
-}
-
-double _CoeffDataAnisotropy::_R2Dsafe(double q, double p)
-{
- if (p <= q && p >= 0 && p < std::numeric_limits<double>::infinity()) {
- if (q == 0) return 0;
- if (q == std::numeric_limits<double>::infinity()) {
- Msg::Warning("q had infinite value in computation of 2D Raniso");
- return 1;
- }
- return (q-p) / (q+p);
- }
- else {
- Msg::Error("Wrong argument for 2d Raniso (%g, %g)", q, p);
- return -1;
- }
-}
-
-double _CoeffDataAnisotropy::_R2Dsafe(double a)
-{
- if (a >= 1) {
- if (a == std::numeric_limits<double>::infinity())
- return 1;
- return (a - 1) / (a + 1);
- }
- else {
- Msg::Error("Wrong argument for 2d metRanisoric (%g)", a);
- return -1;
- }
-}
-
-double _CoeffDataAnisotropy::_R3Dsafe(double q, double p, double J)
-{
- if (p < 0 || p == std::numeric_limits<double>::infinity() || J != J) {
- Msg::Error("Wrong argument for 3d Raniso (%g, %g, %g)", q, p, J);
- return -1;
- }
-
- if (q == 0) return 0;
- if (q > 0 && p == 0) return 1;
- if (q == std::numeric_limits<double>::infinity()) {
- Msg::Warning("q had infinite value in computation of 3D Raniso");
- return 1;
- }
-
- const double a = q/p;
- const double K = J*J/p/p/p;
- return _R3Dsafe(a, K);
-}
-
-double _CoeffDataAnisotropy::_R3Dsafe(double a, double K)
-{
- if (a < 1 || K < 0) {
- if (K >= 0 && a > 1-1e-6 && K < 1e-12) {
- return 0;
- }
- Msg::Error("Wrong argument for 3d Raniso (%g, %g)", a, K);
- _CoeffDataAnisotropy::youReceivedAnError();
- return -1;
- }
-
- if (a == std::numeric_limits<double>::infinity() ||
- K == std::numeric_limits<double>::infinity()) {
- Msg::Warning("a and/or K had infinite value in computation of 3D Raniso");
- return 1;
- }
-
- const double w = _wSafe(a, K);
-
- const double phi = std::acos(w) / 3;
- const double R = (a + 2*std::cos(phi + 2*M_PI/3)) / (a + 2*std::cos(phi));
- return std::max(.0, std::min(1., R));
-}
-
-double _CoeffDataAnisotropy::_wSafe(double a, double K) {
- if (a >= 1 && K >= 0) {
- if (a == std::numeric_limits<double>::infinity() ||
- K == std::numeric_limits<double>::infinity()) {
- Msg::Error("Cannot compute w with infinite numbers (%g, %g)", a, K);
- return 0;
- }
- const double w = _w(a, K);
- if (w > 1) {
- if (w < 1+1e-5 || w < 1+K*1e-14) return 1;
- else {
- Msg::Error("outside the domain w(%g, %g) = %g", a, K, w);
- hasError = true;
- }
- }
- else if (w < -1) {
- if (w > -1-1e-5 || w > -1-K*1e-14) return -1;
- else {
- Msg::Error("outside the domain w(%g, %g) = %g", a, K, w);
- hasError = true;
- }
- }
- return w;
- }
- else {
- youReceivedAnError();
- Msg::Error("Wrong argument for w (%g, %g)", a, K);
- return -1;
- }
-}
-
-int _CoeffDataAnisotropy::metricOrder(MElement *el)
-{
- const int order = el->getPolynomialOrder();
-
- switch (el->getType()) {
- case TYPE_PNT : return 0;
-
- case TYPE_LIN : return order;
-
- case TYPE_TRI :
- case TYPE_TET :
- case TYPE_PYR : return 2*order-2;
-
- case TYPE_QUA :
- case TYPE_PRI :
- case TYPE_HEX : return 2*order;
-
- default :
- Msg::Error("Unknown element type %d, returning -1", el->getType());
- return -1;
- }
-}
-
-void _CoeffDataAnisotropy::interpolateForMatlab(const fullMatrix<double> &nodes) const
-{
- fullVector<double> jac;
- fullMatrix<double> metric;
- _bfsMet->interpolate(_coeffsMetric, nodes, metric);
-
- switch (_coeffsMetric.size2()) {
- case 1:
- // nothing to do
- break;
-
- case 3:
- _CoeffDataAnisotropy::file << "q p myR eigR" << std::endl;
- _CoeffDataAnisotropy::file << _bfsMet->getNumLagCoeff();
- _CoeffDataAnisotropy::file << " " << metric.size1() << std::endl;
- for (int i = 0; i < metric.size1(); ++i) {
- // Compute from q, p
- double p = pow(metric(i, 1), 2);
- p += pow(metric(i, 2), 2);
- p = std::sqrt(p);
- double myR = std::sqrt(_R2Dsafe(metric(i, 0), p));
- // Comppute from tensor
- fullMatrix<double> metricTensor(2, 2);
- metricTensor(0, 0) = metric(i, 0) + metric(i, 1);
- metricTensor(1, 1) = metric(i, 0) - metric(i, 1);
- metricTensor(0, 1) = metricTensor(1, 0) = metric(i, 2);
- fullVector<double> valReal(2), valImag(2);
- fullMatrix<double> vecLeft(2, 2), vecRight(2, 2);
- metricTensor.eig(valReal, valImag, vecLeft, vecRight, true);
- double eigR = std::sqrt(valReal(0) / valReal(1));
- _CoeffDataAnisotropy::file << metric(i, 0) << " ";
- _CoeffDataAnisotropy::file << p << " ";
- _CoeffDataAnisotropy::file << myR << " ";
- _CoeffDataAnisotropy::file << eigR << std::endl;
- }
- break;
-
- case 7:
- _CoeffDataAnisotropy::file << "a K myR eigR" << std::endl;
- _CoeffDataAnisotropy::file << _bfsMet->getNumLagCoeff();
- _CoeffDataAnisotropy::file << " " << metric.size1() << std::endl;
- _bfsDet->interpolate(_coeffsJacDet, nodes, jac);
- for (int i = 0; i < metric.size1(); ++i) {
- // Compute from q, p, J
- double p = 0;
- for (int k = 1; k < 7; ++k) p += pow(metric(i, k), 2);
- p = std::sqrt(p);
- double myR = std::sqrt(_R3Dsafe(metric(i, 0), p, jac(i)));
- // Compute from tensor
- fullMatrix<double> metricTensor(3, 3);
- for (int k = 0; k < 3; ++k) {
- static double fact1 = std::sqrt(6.);
- static double fact2 = std::sqrt(3.);
- const int ki = k%2;
- const int kj = std::min(k+1, 2);
- metricTensor(k, k) = metric(i, k+1)*fact1 + metric(i, 0);
- metricTensor(ki, kj) = metricTensor(kj, ki) = metric(i, k+4)*fact2;
- }
- fullVector<double> valReal(3), valImag(3);
- fullMatrix<double> vecLeft(3, 3), vecRight(3, 3);
- metricTensor.eig(valReal, valImag, vecLeft, vecRight, true);
- double eigR = std::sqrt(valReal(0) / valReal(2));
- _CoeffDataAnisotropy::file << metric(i, 0)/p << " ";
- _CoeffDataAnisotropy::file << jac(i)*jac(i)/p/p/p << " ";
- _CoeffDataAnisotropy::file << myR << " ";
- _CoeffDataAnisotropy::file << eigR << std::endl;
- }
- break;
-
- default:
- Msg::Error("Wrong number of functions for metric: %d", _coeffsMetric.size2());
- }
-}
-
-void _CoeffDataAnisotropy::statsForMatlab(MElement *el, int num) const
-{
- if (num > 1000000) return;
- static unsigned int aa = 0;
- if (++aa < 1000) {
- std::stringstream name;
- name << "metricStat_" << el->getNum() << "_";
- name << (num % 10);
- name << (num % 100)/10;
- name << (num % 1000)/100;
- name << (num % 10000)/1000;
- name << (num % 100000)/10000;
- name << ".txt";
- _CoeffDataAnisotropy::file.open(name.str().c_str(), std::fstream::out);
-
- {
- _CoeffDataAnisotropy::file << "mina minKfast minKsharp cmin beta_m beta_M beta c_m a_int K_int pdRda pdRdK a0" << std::endl;
- double mina, minKs, minKf, gamma, beta_m, beta_M, beta, c_m;
- double a_int, K_int, p_dRda, p_dRdK, a0;
-
- mina = _computeLowerBoundA();
-
- if (el->getDim() == 3) {
- minKs = minKf = _computeLowerBoundK();
- double minK = minKf;
- _moveInsideDomain(mina, minK, true);
- _computePseudoDerivatives(mina, minK, p_dRda, p_dRdK);
- double slope = -p_dRda/p_dRdK;
- gamma = .5*(3*minK/mina - slope);
- _computeBoundingCurve(beta_m, gamma, true);
- _computeBoundingCurve(beta_M, gamma, false);
- beta = -p_dRda/(p_dRdK*3*mina*mina);
- c_m = -1;
- if (p_dRda < 0) {
- a_int = mina, K_int = minK;
- _intersectionCurveLeftCorner(beta_m, gamma, a_int, K_int);
- }
- else {
- a_int = mina, K_int = minK;
- _intersectionCurveLeftCorner(beta_M, gamma, a_int, K_int);
- }
- a0 = _computeAbscissaMinR(mina, minK);
- }
- else {
- minKs = minKf = gamma = beta_m = beta_M = beta = c_m = -1;
- a_int = K_int = p_dRda = p_dRdK = a0 = -1;
- }
-
- _CoeffDataAnisotropy::file << std::setprecision(15);
- _CoeffDataAnisotropy::file << mina << " " << minKf << " " << minKs << " " << std::endl;
- _CoeffDataAnisotropy::file << gamma << " " << beta_m << " " << beta_M << " " << std::endl;
- _CoeffDataAnisotropy::file << beta << " " << c_m << std::endl;
- _CoeffDataAnisotropy::file << a_int << " " << K_int << std::endl;
- _CoeffDataAnisotropy::file << p_dRda << " " << p_dRdK << std::endl;
- _CoeffDataAnisotropy::file << a0 << std::endl;
- _CoeffDataAnisotropy::file << std::setprecision(8);
- }
- }
-
- fullMatrix<double> samplingPoints;
- bool serendip = false;
- gmshGeneratePoints(FuncSpaceData(el, 20, &serendip), samplingPoints);
- interpolateForMatlab(samplingPoints);
- _CoeffDataAnisotropy::file.close();
-}
-
// Miscellaneous
template<typename Comp>
void _subdivideDomainsMinOrMax(std::vector<_CoeffData*> &domains,
@@ -1890,8 +976,6 @@ void _subdivideDomainsMinOrMax(std::vector<_CoeffData*> &domains,
}
}
if (k > 1000) {
- if (domains[0]->getNumMeasure() == 3)
- _CoeffDataAnisotropy::youReceivedAnError();
Msg::Error("Max subdivision (1000) (size %d)", domains.size());
}
}
@@ -1921,7 +1005,6 @@ double _computeBoundRational(const fullVector<double> &numerator,
Msg::Fatal("In order to compute a bound on a rational function, I need "
"vectors of the same size! (%d vs %d)", numerator.size(),
denominator.size());
- _CoeffDataAnisotropy::youReceivedAnError();
return 0;
}
diff --git a/Mesh/qualityMeasuresJacobian.h b/Mesh/qualityMeasuresJacobian.h
index fcc7d52a92e3274b8d7ec60246aaa9183e8736a3..2ed23400530756962ad0bc52b2c3ff5d2875acc7 100644
--- a/Mesh/qualityMeasuresJacobian.h
+++ b/Mesh/qualityMeasuresJacobian.h
@@ -20,10 +20,6 @@ void minMaxJacobianDeterminant(MElement *el, double &min, double &max,
double minScaledJacobian(MElement *el,
bool knownValid = false,
bool reversedOk = false);
-double minAnisotropyMeasure(MElement *el,
- bool knownValid = false,
- bool reversedOk = false,
- int n = 5); //n is fordebug
double minIsotropyMeasure(MElement *el,
bool knownValid = false,
bool reversedOk = false);
@@ -110,79 +106,6 @@ private:
fullMatrix<double> &coeffScaledJacobian) const;
};
-class _CoeffDataAnisotropy: public _CoeffData
-{
-private:
- const fullVector<double> _coeffsJacDet;
- const fullMatrix<double> _coeffsMetric;
- const bezierBasis *_bfsDet, *_bfsMet;
- MElement *_elForDebug;
- int _numForDebug;
- static bool hasError;//fordebug
- static std::fstream file;//fordebug
-
-public:
- static std::vector<double> mytimes;//fordebug
- static std::vector<int> mynumbers;//fordebug
- static MElement *currentElement;//fordebug
- static _CoeffDataAnisotropy *current;//fordebug
- _CoeffDataAnisotropy(fullVector<double> &det,
- fullMatrix<double> &metric,
- const bezierBasis *bfsDet,
- const bezierBasis *bfsMet,
- int depth, MElement *,
- int num = 0);
- ~_CoeffDataAnisotropy() {}
-
- bool boundsOk(double minL, double maxL) const;
- void getSubCoeff(std::vector<_CoeffData*>&) const;
- int getNumMeasure() const {return 3;}//fordebug
-
- static int metricOrder(MElement *el);
- static void fillMetricCoeff(const GradientBasis*,
- const fullMatrix<double> &nodes,
- fullMatrix<double> &coeff,
- int dim, bool ideal);
-
- static bool noMoreComputationPlease() {return hasError;}//fordebug
- static void youReceivedAnError()
- {
- current->statsForMatlab(currentElement, 20);
- hasError = true;
- }//fordebug
-
-private:
- void _computeAtCorner(double &min, double &max) const;
- double _computeLowerBound() const;
-
- double _computeLowerBoundA() const;
- double _computeLowerBoundK() const;
- void _computeBoundingCurve(double &beta, double gamma, bool lowerBound) const;
-
- static bool _moveInsideDomain(double &a, double &K, bool bottomleftCorner);
- static bool _computePseudoDerivatives(double a, double K,
- double &dRda, double &dRdK);
- static bool _computeDerivatives(double a, double K,
- double &dRda, double &dRdK,
- double &dRdaa, double &dRdaK, double &dRdKK);
- static int _intersectionCurveLeftCorner(double beta, double gamma,
- double &mina, double &minK);
- static double _computeAbscissaMinR(double aStart, double K);
-
- static double _R2Dsafe(double q, double p);
- static double _R2Dsafe(double a);
- static double _R3Dsafe(double q, double p, double J);
- static double _R3Dsafe(double a, double K);
- static double _wSafe(double a, double K);
- static inline double _w(double a, double K) {
- return .5 * (K + (3-a*a)*a);
- }
-
-public:
- void interpolateForMatlab(const fullMatrix<double> &nodes) const;//fordebug
- void statsForMatlab(MElement *el, int) const;//fordebug
-};
-
class _CoeffDataIsotropy: public _CoeffData
{
private:
diff --git a/Mesh/yamakawa.cpp b/Mesh/yamakawa.cpp
index 76cfd27ca9fed573c88ca50e7a910faff722ad04..a4f852b3ef4e84d7e122f721eb3c11cd155dc256 100644
--- a/Mesh/yamakawa.cpp
+++ b/Mesh/yamakawa.cpp
@@ -5611,12 +5611,12 @@ void PostOp::split_pyramids(GRegion* gr){
tempC3->getVolumeSign(),
tempC4->getVolumeSign());
}
- double qA = tempA1->minAnisotropyMeasure() + tempA2->minAnisotropyMeasure() +
- tempA3->minAnisotropyMeasure() + tempA4->minAnisotropyMeasure();
- double qB = tempB1->minAnisotropyMeasure() + tempB2->minAnisotropyMeasure() +
- tempB3->minAnisotropyMeasure() + tempB4->minAnisotropyMeasure();
- double qC = tempC1->minAnisotropyMeasure() + tempC2->minAnisotropyMeasure() +
- tempC3->minAnisotropyMeasure() + tempC4->minAnisotropyMeasure();
+ double qA = tempA1->minIsotropyMeasure() + tempA2->minIsotropyMeasure() +
+ tempA3->minIsotropyMeasure() + tempA4->minIsotropyMeasure();
+ double qB = tempB1->minIsotropyMeasure() + tempB2->minIsotropyMeasure() +
+ tempB3->minIsotropyMeasure() + tempB4->minIsotropyMeasure();
+ double qC = tempC1->minIsotropyMeasure() + tempC2->minIsotropyMeasure() +
+ tempC3->minIsotropyMeasure() + tempC4->minIsotropyMeasure();
if (qA > qB && qA > qC) {
Msg::Info("A");
gr->addTetrahedron(tempA1);