diff --git a/Mesh/qualityMeasuresJacobian.cpp b/Mesh/qualityMeasuresJacobian.cpp index d48d3206393c11771d9d9cb73de9eee98d131bab..f70d472e6028f4e0ad2b84188a6c752da8a3398a 100644 --- a/Mesh/qualityMeasuresJacobian.cpp +++ b/Mesh/qualityMeasuresJacobian.cpp @@ -663,10 +663,18 @@ void _CoeffDataScaledJac::_computeAtCorner(double &min, double &max) const 1/v(i,1)/v(i,2)) / 3; break; case TYPE_TET: - sJ = cTet * _coeffsJacDet(i) * (1/v(i,0)/v(i,1)/v(i,2) + - 1/v(i,0)/v(i,3)/v(i,4) + - 1/v(i,1)/v(i,3)/v(i,5) + - 1/v(i,2)/v(i,4)/v(i,5)) / 4; + sJ = cTet * _coeffsJacDet(i) * (1/v(i,0)/v(i,5)/v(i,1) + + 1/v(i,0)/v(i,5)/v(i,2) + + 1/v(i,0)/v(i,5)/v(i,3) + + 1/v(i,0)/v(i,5)/v(i,4) + + 1/v(i,1)/v(i,4)/v(i,0) + + 1/v(i,1)/v(i,4)/v(i,2) + + 1/v(i,1)/v(i,4)/v(i,3) + + 1/v(i,1)/v(i,4)/v(i,5) + + 1/v(i,2)/v(i,3)/v(i,0) + + 1/v(i,2)/v(i,3)/v(i,1) + + 1/v(i,2)/v(i,3)/v(i,4) + + 1/v(i,2)/v(i,3)/v(i,5)) / 12; break; case TYPE_PRI: sJ = cTri * _coeffsJacDet(i) * (1/v(i,0)/v(i,1)/v(i,2) + @@ -743,15 +751,39 @@ double _CoeffDataScaledJac::_computeLowerBound() const return cTri*result/3; case TYPE_TET: - raiser->computeCoeff(prox[0], prox[1], prox[2], coeffDenominator); - result += _computeBoundRational(_coeffsJacDet, coeffDenominator, true); - raiser->computeCoeff(prox[0], prox[3], prox[4], coeffDenominator); - result += _computeBoundRational(_coeffsJacDet, coeffDenominator, true); - raiser->computeCoeff(prox[1], prox[3], prox[5], coeffDenominator); - result += _computeBoundRational(_coeffsJacDet, coeffDenominator, true); - raiser->computeCoeff(prox[2], prox[4], prox[5], coeffDenominator); - result += _computeBoundRational(_coeffsJacDet, coeffDenominator, true); - return cTet*result/4; + { + fullVector<double> thirdTerm, coeffNum1, tmp; + thirdTerm = prox[1]; + thirdTerm.axpy(prox[2]); + thirdTerm.axpy(prox[3]); + thirdTerm.axpy(prox[4]); + raiser->computeCoeff(prox[0], prox[5], thirdTerm, coeffNum1); + thirdTerm = prox[0]; + thirdTerm.axpy(prox[2]); + thirdTerm.axpy(prox[3]); + thirdTerm.axpy(prox[5]); + raiser->computeCoeff(prox[1], prox[4], thirdTerm, tmp); + coeffNum1.axpy(tmp); + thirdTerm = prox[0]; + thirdTerm.axpy(prox[1]); + thirdTerm.axpy(prox[4]); + thirdTerm.axpy(prox[5]); + raiser->computeCoeff(prox[2], prox[3], thirdTerm, tmp); + coeffNum1.axpy(tmp); + + fullVector<double> coeffDen1, coeffDen2; + raiser->computeCoeff(prox[0], prox[1], prox[2], coeffDen1); + raiser->computeCoeff(prox[3], prox[4], prox[5], coeffDen2); + + fullVector<double> &coeffNumerator = tmp; + bezierBasisRaiser *raiserBis; + raiserBis = raiser->getRaisedBezierBasis(3)->getRaiser(); + raiserBis->computeCoeff(coeffNum1, _coeffsJacDet, coeffNumerator); + raiserBis->computeCoeff(coeffDen1, coeffDen2, coeffDenominator); + + result = _computeBoundRational(coeffNumerator, coeffDenominator, true); + return cTet*result/12; + } case TYPE_PYR: raiser->computeCoeff(prox[0], prox[1], prox[2], coeffDenominator); diff --git a/Numeric/FuncSpaceData.h b/Numeric/FuncSpaceData.h index b1bb7e57bd2c629ad96cc7718fbcbcbc3fb7636d..637c707428a63b4894adceb72035522300c36ff2 100644 --- a/Numeric/FuncSpaceData.h +++ b/Numeric/FuncSpaceData.h @@ -42,7 +42,7 @@ public: : _tag(-1), _spaceOrder(-1), _serendipity(false), _nij(-1), _nk(-1), _pyramidalSpace(false) {} - // Constructors of the function space of a different order + // Constructors for the function space of a different order FuncSpaceData(const FuncSpaceData &fsd, int order, const bool *serendip = NULL); diff --git a/Numeric/bezierBasis.cpp b/Numeric/bezierBasis.cpp index 88c5339af30ce46564fa6a54e1ccfc989b55179c..805f6e2154f252318a2273428b848f6e4771b16e 100644 --- a/Numeric/bezierBasis.cpp +++ b/Numeric/bezierBasis.cpp @@ -959,6 +959,19 @@ void bezierBasisRaiser::_fillRaiserDataPyr() } } +const bezierBasis* bezierBasisRaiser::getRaisedBezierBasis(int mult2or3) +{ + if (mult2or3 != 2 && mult2or3 != 3) { + Msg::Error("There is no reason that I give you the bezier basis of order " + "%d times greater. Ask to another class/object", mult2or3); + return NULL; + } + FuncSpaceData initFuncSpace = _bfs->getFuncSpaceData(); + int newOrder = initFuncSpace.spaceOrder() * mult2or3; + FuncSpaceData raisedFuncSpace(initFuncSpace, newOrder); + return BasisFactory::getBezierBasis(raisedFuncSpace); +} + void bezierBasisRaiser::computeCoeff(const fullVector<double> &coeffA, const fullVector<double> &coeffB, fullVector<double> &coeffSquare) diff --git a/Numeric/bezierBasis.h b/Numeric/bezierBasis.h index c0d50bc70754fbfaa198df39ba5dfff574c41f70..ac31eb01a0c0a184645a36f1b7dcd6a1aa540c60 100644 --- a/Numeric/bezierBasis.h +++ b/Numeric/bezierBasis.h @@ -102,6 +102,8 @@ public: _fillRaiserData(); }; + const bezierBasis* getRaisedBezierBasis(int multipliedBy2or3); + // const bezierBasis* getRaisedBezierBasis(int raised) const; void computeCoeff(const fullVector<double> &coeffA,