From af47f665ec73dc9e5c627dc335fcb2d946b12835 Mon Sep 17 00:00:00 2001
From: Amaury Johnan <amjohnen@gmail.com>
Date: Wed, 21 Sep 2016 14:15:16 +0000
Subject: [PATCH] update scaled jacobian for tets
---
Mesh/qualityMeasuresJacobian.cpp | 58 +++++++++++++++++++++++++-------
Numeric/FuncSpaceData.h | 2 +-
Numeric/bezierBasis.cpp | 13 +++++++
Numeric/bezierBasis.h | 2 ++
4 files changed, 61 insertions(+), 14 deletions(-)
diff --git a/Mesh/qualityMeasuresJacobian.cpp b/Mesh/qualityMeasuresJacobian.cpp
index d48d320639..f70d472e60 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 b1bb7e57bd..637c707428 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 88c5339af3..805f6e2154 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 c0d50bc707..ac31eb01a0 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,
--
GitLab