diff --git a/Mesh/qualityMeasuresJacobian.cpp b/Mesh/qualityMeasuresJacobian.cpp index 8dfc942ffa040d04e9c2208b6a2c22d759d8325e..1fedcb7d096284188935c4941d62a73e8c447a15 100644 --- a/Mesh/qualityMeasuresJacobian.cpp +++ b/Mesh/qualityMeasuresJacobian.cpp @@ -35,7 +35,7 @@ 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); +double _CoeffDataScaledJac::cPyr = std::sqrt(2)*4; void minMaxJacobianDeterminant(MElement *el, double &min, double &max) { @@ -115,8 +115,8 @@ double minScaledJacobian(MElement *el, bool knownValid, bool reversedOk) jacDetSpace = FuncSpaceData(el, jacOrder, &serendipFalse); break; case TYPE_PYR: - jacMatSpace = FuncSpaceData(el, false, order, order, &serendipFalse); - jacDetSpace = FuncSpaceData(el, false, jacOrder, jacOrder, &serendipFalse); + jacMatSpace = FuncSpaceData(el, false, order, order-1, &serendipFalse); + jacDetSpace = FuncSpaceData(el, false, jacOrder, jacOrder-3, &serendipFalse); break; default: Msg::Error("Scaled Jacobian not implemented for type of element %d", el->getType()); @@ -154,6 +154,12 @@ double minScaledJacobian(MElement *el, bool knownValid, bool reversedOk) _subdivideDomains(domains); +// Msg::Info("element %d, type %d, numSub %d", el->getNum(), el->getType(), +// domains.size()/7); + if (domains.size()/7 > 500) {//fordebug + Msg::Info("S 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) { @@ -249,8 +255,8 @@ double minAnisotropyMeasure(MElement *el, ); _subdivideDomains(domains); - if (domains.size()/7+1 > 10) {//fordebug - Msg::Info("too much subdivision: %d (el %d)", domains.size()/7+1, el->getNum()); + 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(); @@ -267,6 +273,104 @@ double minAnisotropyMeasure(MElement *el, return min; } +double minIsotropyMeasure(MElement *el, + bool knownValid, + bool reversedOk) +{ + bool isReversed = false; + if (!knownValid) { + double jmin, jmax; + minMaxJacobianDeterminant(el, jmin, jmax); + if (jmax < 0) { + if (!reversedOk) return 0; + isReversed = true; + } + else if (jmin <= 0) return 0; + } + // Computation of the minimum (or maximum) of the scaled Jacobian should never + // be performed to invalid elements (for which the measure is 0). + + fullMatrix<double> nodesXYZ(el->getNumVertices(), 3); + el->getNodesCoord(nodesXYZ); + + const JacobianBasis *jacBasis; + const GradientBasis *gradBasis; + + const int type = el->getType(); + const int order = el->getPolynomialOrder(); + const int jacOrder = order * el->getDim(); + const bool serendipFalse = false; + + FuncSpaceData jacMatSpace, jacDetSpace; + + switch(type) { + case TYPE_TRI: + jacMatSpace = FuncSpaceData(el, order-1, &serendipFalse); + jacDetSpace = FuncSpaceData(el, jacOrder-2, &serendipFalse); + break; + case TYPE_TET: + jacMatSpace = FuncSpaceData(el, order-1, &serendipFalse); + jacDetSpace = FuncSpaceData(el, jacOrder-3, &serendipFalse); + break; + case TYPE_QUA: + case TYPE_HEX: + case TYPE_PRI: + jacMatSpace = FuncSpaceData(el, order, &serendipFalse); + jacDetSpace = FuncSpaceData(el, jacOrder, &serendipFalse); + break; + case TYPE_PYR: + jacMatSpace = FuncSpaceData(el, false, order, order-1, &serendipFalse); + jacDetSpace = FuncSpaceData(el, false, jacOrder, jacOrder-3, &serendipFalse); + break; + default: + Msg::Error("Scaled Jacobian not implemented for type of element %d", el->getType()); + return -1; + } + gradBasis = BasisFactory::getGradientBasis(jacMatSpace); + jacBasis = BasisFactory::getJacobianBasis(jacDetSpace); + + fullVector<double> coeffDetBez; + { + fullVector<double> coeffDetLag(jacBasis->getNumJacNodes()); + jacBasis->getSignedIdealJacobian(nodesXYZ, coeffDetLag); + + coeffDetBez.resize(jacBasis->getNumJacNodes()); + jacBasis->lag2Bez(coeffDetLag, coeffDetBez); + + if (isReversed) coeffDetBez.scale(-1); + } + + fullMatrix<double> coeffMatBez; + { + fullMatrix<double> coeffMatLag(gradBasis->getNumSamplingPoints(), 9); + gradBasis->getAllIdealGradientsFromNodes(nodesXYZ, coeffMatLag); + + coeffMatBez.resize(gradBasis->getNumSamplingPoints(), 9); + gradBasis->lag2Bez(coeffMatLag, coeffMatBez); + if (el->getDim() == 2) coeffMatBez.resize(coeffMatBez.size1(), 6, false); + } + + std::vector<_CoeffData*> domains; + domains.push_back( + new _CoeffDataIsotropy(coeffDetBez, coeffMatBez, jacBasis->getBezier(), + gradBasis->getBezier(), 0) + ); + + _subdivideDomains(domains); + if (domains.size()/7 > 500) {//fordebug + Msg::Info("I 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]; + } + return min; +} + //double minSampledAnisotropyMeasure(MElement *el, int deg, bool writeInFile)//fordebug //{ // fullMatrix<double> samplingPoints; @@ -374,9 +478,42 @@ _CoeffDataScaledJac::_CoeffDataScaledJac(fullVector<double> &det, _computeAtCorner(_minL, _maxL); _minB = _computeLowerBound(); +// Msg::Info("%g %g %g", _minB, _minL, _maxL);//fordebug +// _coeffsJacDet.print("_coeffsJacDet"); +// _coeffsJacMat.print("_coeffsJacMat"); // computation of _maxB not implemented for now } +bool _CoeffDataScaledJac::boundsOk(double minL, double maxL) const +{ + static double tol = 1e-3; + return minL - _minB < tol; +} + +void _CoeffDataScaledJac::getSubCoeff(std::vector<_CoeffData*> &v) const +{ + v.clear(); + v.reserve(_bfsDet->getNumDivision()); + fullVector<double> subCoeffD; + fullMatrix<double> subCoeffM; + _bfsDet->subdivideBezCoeff(_coeffsJacDet, subCoeffD); + _bfsMat->subdivideBezCoeff(_coeffsJacMat, subCoeffM); + + int szD = _coeffsJacDet.size(); + int szM1 = _coeffsJacMat.size1(); + int szM2 = _coeffsJacMat.size2(); + for (int i = 0; i < _bfsDet->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); + _CoeffDataScaledJac *newData; + newData = new _CoeffDataScaledJac(coeffD, coeffM, _bfsDet, _bfsMat, + _depth+1, _type); + v.push_back(newData); + } +} + void _CoeffDataScaledJac::_computeAtCorner(double &min, double &max) const { min = std::numeric_limits<double>::infinity(); @@ -535,79 +672,226 @@ void _CoeffDataScaledJac::_getCoeffLengthVectors(fullMatrix<double> &coeff, const fullMatrix<double> &mat = _coeffsJacMat; - for (int i = 0; i < sz1; i++) { - coeff(i, 0) = std::sqrt(pow_int(mat(i, 0), 2) + - pow_int(mat(i, 1), 2) + - pow_int(mat(i, 2), 2) ); - coeff(i, 1) = std::sqrt(pow_int(mat(i, 3), 2) + - pow_int(mat(i, 4), 2) + - pow_int(mat(i, 5), 2) ); - + if (_type != TYPE_PYR) { + for (int i = 0; i < sz1; i++) { + coeff(i, 0) = std::sqrt(pow_int(mat(i, 0), 2) + + pow_int(mat(i, 1), 2) + + pow_int(mat(i, 2), 2) ); + coeff(i, 1) = std::sqrt(pow_int(mat(i, 3), 2) + + pow_int(mat(i, 4), 2) + + pow_int(mat(i, 5), 2) ); + } if (mat.size2() > 6) { - coeff(i, 2) = std::sqrt(pow_int(mat(i, 6), 2) + - pow_int(mat(i, 7), 2) + - pow_int(mat(i, 8), 2) ); + for (int i = 0; i < sz1; i++) { + coeff(i, 2) = std::sqrt(pow_int(mat(i, 6), 2) + + pow_int(mat(i, 7), 2) + + pow_int(mat(i, 8), 2) ); + } } else if (_type == TYPE_TRI) { - coeff(i, 2) = std::sqrt(pow_int(mat(i, 3) - mat(i, 0), 2) + - pow_int(mat(i, 4) - mat(i, 1), 2) + - pow_int(mat(i, 5) - mat(i, 2), 2) ); + for (int i = 0; i < sz1; i++) { + coeff(i, 2) = std::sqrt(pow_int(mat(i, 3) - mat(i, 0), 2) + + pow_int(mat(i, 4) - mat(i, 1), 2) + + pow_int(mat(i, 5) - mat(i, 2), 2) ); + } } - switch (_type) { case TYPE_TET: case TYPE_PRI: - coeff(i, 3) = std::sqrt(pow_int(mat(i, 3) - mat(i, 0), 2) + - pow_int(mat(i, 4) - mat(i, 1), 2) + - pow_int(mat(i, 5) - mat(i, 2), 2) ); - break; - case TYPE_PYR: - coeff(i, 3) = std::sqrt(pow_int(mat(i, 6) - mat(i, 0) - mat(i, 3), 2) + + for (int i = 0; i < sz1; i++) { + coeff(i, 3) = std::sqrt(pow_int(mat(i, 3) - mat(i, 0), 2) + + pow_int(mat(i, 4) - mat(i, 1), 2) + + pow_int(mat(i, 5) - mat(i, 2), 2) ); + } + } + if (_type == TYPE_TET) { + for (int i = 0; i < sz1; i++) { + coeff(i, 4) = std::sqrt(pow_int(mat(i, 6) - mat(i, 0), 2) + + pow_int(mat(i, 7) - mat(i, 1), 2) + + pow_int(mat(i, 8) - mat(i, 2), 2) ); + coeff(i, 5) = std::sqrt(pow_int(mat(i, 6) - mat(i, 3), 2) + + pow_int(mat(i, 7) - mat(i, 4), 2) + + pow_int(mat(i, 8) - mat(i, 5), 2) ); + } + } + } + else { + for (int i = 0; i < sz1; i++) { + coeff(i, 0) = std::sqrt(pow_int(2*mat(i, 0), 2) + + pow_int(2*mat(i, 1), 2) + + pow_int(2*mat(i, 2), 2) ); + coeff(i, 1) = std::sqrt(pow_int(2*mat(i, 3), 2) + + pow_int(2*mat(i, 4), 2) + + pow_int(2*mat(i, 5), 2) ); + coeff(i, 2) = std::sqrt(pow_int(mat(i, 6) + mat(i, 0) + mat(i, 3), 2) + + pow_int(mat(i, 7) + mat(i, 1) + mat(i, 4), 2) + + pow_int(mat(i, 8) + mat(i, 2) + mat(i, 5), 2) ); + coeff(i, 3) = std::sqrt(pow_int(mat(i, 6) - mat(i, 0) + mat(i, 3), 2) + + pow_int(mat(i, 7) - mat(i, 1) + mat(i, 4), 2) + + pow_int(mat(i, 8) - mat(i, 2) + mat(i, 5), 2) ); + coeff(i, 4) = std::sqrt(pow_int(mat(i, 6) - mat(i, 0) - mat(i, 3), 2) + pow_int(mat(i, 7) - mat(i, 1) - mat(i, 4), 2) + pow_int(mat(i, 8) - mat(i, 2) - mat(i, 5), 2) ); + coeff(i, 5) = std::sqrt(pow_int(mat(i, 6) + mat(i, 0) - mat(i, 3), 2) + + pow_int(mat(i, 7) + mat(i, 1) - mat(i, 4), 2) + + pow_int(mat(i, 8) + mat(i, 2) - mat(i, 5), 2) ); } + } +} - if (_type == TYPE_TET || _type == TYPE_PYR) { - coeff(i, 4) = std::sqrt(pow_int(mat(i, 6) - mat(i, 0), 2) + - pow_int(mat(i, 7) - mat(i, 1), 2) + - pow_int(mat(i, 8) - mat(i, 2), 2) ); - coeff(i, 5) = std::sqrt(pow_int(mat(i, 6) - mat(i, 3), 2) + - pow_int(mat(i, 7) - mat(i, 4), 2) + - pow_int(mat(i, 8) - mat(i, 5), 2) ); - } +// Isotropy measure +_CoeffDataIsotropy::_CoeffDataIsotropy(fullVector<double> &det, + fullMatrix<double> &mat, + const bezierBasis *bfsDet, + const bezierBasis *bfsMat, + int depth) +: _CoeffData(depth), _coeffsJacDet(det.getDataPtr(), det.size()), + _coeffsJacMat(mat.getDataPtr(), mat.size1(), mat.size2()), + _bfsDet(bfsDet), _bfsMat(bfsMat) +{ + if (!det.getOwnData() || !mat.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); + mat.setOwnData(false); + const_cast<fullVector<double>&>(_coeffsJacDet).setOwnData(true); + const_cast<fullMatrix<double>&>(_coeffsJacMat).setOwnData(true); + + _computeAtCorner(_minL, _maxL); + + _minB = 0; + if (boundsOk(_minL, _maxL)) return; + else _minB = _computeLowerBound(); +// Msg::Info("%g %g %g ", _minB, _minL, _maxL); +// _coeffsJacDet.print("_coeffsJacDet"); +// _coeffsJacMat.print("_coeffsJacMat"); + + // _maxB not used for now } -bool _CoeffDataScaledJac::boundsOk(double minL, double maxL) const +bool _CoeffDataIsotropy::boundsOk(double minL, double maxL) const { static double tol = 1e-3; - return minL - _minB < tol; + return minL < tol*1e-3 || (_minB > minL-tol && _minB > minL*(1-100*tol)); } -void _CoeffDataScaledJac::getSubCoeff(std::vector<_CoeffData*> &v) const +void _CoeffDataIsotropy::getSubCoeff(std::vector<_CoeffData*> &v) const { v.clear(); - v.reserve(_bfsDet->getNumDivision()); - fullVector<double> subCoeffD; + v.reserve(_bfsMat->getNumDivision()); fullMatrix<double> subCoeffM; - _bfsDet->subdivideBezCoeff(_coeffsJacDet, subCoeffD); + fullVector<double> subCoeffD; _bfsMat->subdivideBezCoeff(_coeffsJacMat, subCoeffM); + _bfsDet->subdivideBezCoeff(_coeffsJacDet, subCoeffD); int szD = _coeffsJacDet.size(); int szM1 = _coeffsJacMat.size1(); int szM2 = _coeffsJacMat.size2(); - for (int i = 0; i < _bfsDet->getNumDivision(); i++) { + for (int i = 0; i < _bfsMat->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); - _CoeffDataScaledJac *newData; - newData = new _CoeffDataScaledJac(coeffD, coeffM, _bfsDet, _bfsMat, - _depth+1, _type); + _CoeffDataIsotropy *newData + = new _CoeffDataIsotropy(coeffD, coeffM, _bfsDet, _bfsMat, _depth+1); v.push_back(newData); } } +void _CoeffDataIsotropy::_computeAtCorner(double &min, double &max) const +{ + min = std::numeric_limits<double>::infinity(); + max = -min; + + for (int i = 0; i < _bfsMat->getNumLagCoeff(); i++) { + double p = 0; + for (int k = 0; k < _coeffsJacMat.size2(); ++k) { + p += pow_int(_coeffsJacMat(i, k), 2); + } + double qual; + if (_coeffsJacMat.size2() == 6) + qual = 2 * _coeffsJacDet(i) / p; + else + qual = 3 * std::pow(_coeffsJacDet(i), 2/3.) / p; + min = std::min(min, qual); + max = std::max(max, qual); + } +} + +double _CoeffDataIsotropy::_computeLowerBound() const +{ + // Speedup: If one coeff _coeffsJacDet is negative, we would get + // a negative lower bound. For now, returning 0. + for (int i = 0; i < _coeffsJacDet.size(); ++i) { + if (_coeffsJacDet(i) < 0) { + return 0; + } + } + + // 2D element + if (_coeffsJacMat.size2() == 6) { + fullVector<double> coeffDenominator; + { + bezierBasisRaiser *raiser = _bfsMat->getRaiser(); + fullVector<double> prox; + for (int k = 0; k < _coeffsJacMat.size2(); ++k) { + prox.setAsProxy(_coeffsJacMat, k); + fullVector<double> tmp; + raiser->computeCoeff(prox, prox, tmp); + if (k == 0) coeffDenominator.resize(tmp.size()); + coeffDenominator.axpy(tmp, 1); + } + } + return 2*_computeBoundRational(_coeffsJacDet, coeffDenominator, true); + } + + // 3D element NEW + fullVector<double> coeffDenominator; + { + fullVector<double> P(_coeffsJacMat.size1()); + // element of P are automatically set to 0 + for (int i = 0; i < _coeffsJacMat.size1(); ++i) { + for (int k = 0; k < _coeffsJacMat.size2(); ++k) { + P(i) += _coeffsJacMat(i, k) * _coeffsJacMat(i, k); + } + P(i) = std::sqrt(P(i)); + } + _bfsMat->getRaiser()->computeCoeff(P, P, P, coeffDenominator); + } + + return 3*std::pow(_computeBoundRational(_coeffsJacDet, + coeffDenominator, true), 2/3.); + +// // 3D element OLD +// fullVector<double> coeffNumerator; +// _bfsDet->getRaiser()->computeCoeff(_coeffsJacDet, _coeffsJacDet, coeffNumerator); +// +// fullVector<double> coeffDenominator; +// { +// fullVector<double> coeffDenCbrt; +// bezierBasisRaiser *raiser = _bfsMat->getRaiser(); +// +// fullVector<double> prox; +// for (int k = 0; k < _coeffsJacMat.size2(); ++k) { +// prox.setAsProxy(_coeffsJacMat, k); +// fullVector<double> tmp; +// raiser->computeCoeff(prox, prox, tmp); +// if (k == 0) coeffDenCbrt.resize(tmp.size()); +// coeffDenCbrt.axpy(tmp, 1); +// } +// +// bezierBasisRaiser *raiser2 = raiser->getRaisedBezierBasis(2)->getRaiser(); +// raiser2->computeCoeff(coeffDenCbrt, coeffDenCbrt, coeffDenCbrt, coeffDenominator); +// } +// +// return 3*std::pow(_computeBoundRational(coeffNumerator, +// coeffDenominator, true), 1/3.); +} + // Anisotropy measure _CoeffDataAnisotropy::_CoeffDataAnisotropy(fullVector<double> &det, fullMatrix<double> &metric, @@ -1420,7 +1704,7 @@ void _subdivideDomainsMinOrMax(std::vector<_CoeffData*> &domains, std::vector<_CoeffData*> subs; make_heap(domains.begin(), domains.end(), Comp()); int k = 0; - while (!domains[0]->boundsOk(minL, maxL) && k++ < 100) { + while (!domains[0]->boundsOk(minL, maxL) && k++ < 1000) { _CoeffData *cd = domains[0]; pop_heap(domains.begin(), domains.end(), Comp()); domains.pop_back(); @@ -1434,10 +1718,10 @@ void _subdivideDomainsMinOrMax(std::vector<_CoeffData*> &domains, push_heap(domains.begin(), domains.end(), Comp()); } } - if (k > 100) { + if (k > 1000) { if (domains[0]->getNumMeasure() == 3) _CoeffDataAnisotropy::youReceivedAnError(); - Msg::Error("Max subdivision (100) (size %d)", domains.size()); + Msg::Error("Max subdivision (1000) (size %d)", domains.size()); } } diff --git a/Mesh/qualityMeasuresJacobian.h b/Mesh/qualityMeasuresJacobian.h index 72dc98541413ace50ced46747f7fae10013e3aee..5c4921c60a09dff9010ff4320a706fad424d8007 100644 --- a/Mesh/qualityMeasuresJacobian.h +++ b/Mesh/qualityMeasuresJacobian.h @@ -23,6 +23,9 @@ 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); //double minSampledAnisotropyMeasure(MElement *el, int order,//fordebug // bool writeInFile = false); @@ -174,11 +177,29 @@ public: void statsForMatlab(MElement *el, int) const;//fordebug }; -//inline bool isValid(MElement *el) { -// double min, max; -// minMaxJacobianDeterminant(el, min, max); -// return min > 0; -//} +class _CoeffDataIsotropy: public _CoeffData +{ +private: + const fullVector<double> _coeffsJacDet; + const fullMatrix<double> _coeffsJacMat; + const bezierBasis *_bfsDet, *_bfsMat; + +public: + _CoeffDataIsotropy(fullVector<double> &det, + fullMatrix<double> &metric, + const bezierBasis *bfsDet, + const bezierBasis *bfsMet, + int depth); + ~_CoeffDataIsotropy() {} + + bool boundsOk(double minL, double maxL) const; + void getSubCoeff(std::vector<_CoeffData*>&) const; + int getNumMeasure() const {return 4;}//fordebug + +private: + void _computeAtCorner(double &min, double &max) const; + double _computeLowerBound() const; +}; double _computeBoundRational(const fullVector<double> &numerator, const fullVector<double> &denominator,