Skip to content
Snippets Groups Projects
Commit 5f51edba authored by Amaury Johnen's avatar Amaury Johnen
Browse files

(1) add new 'isotropy' measure: equivalent to the old 'anisotropy' measure but... 

(1) add new 'isotropy' measure: equivalent to the old 'anisotropy' measure but easier to compute. (2) imporvements
parent 70194ee5
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Mesh/qualityMeasuresJacobian.cpp
+ 331
− 47
View file @ 5f51edba
......@@ -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());
}
}
......
Mesh/qualityMeasuresJacobian.h
+ 26
− 5
View file @ 5f51edba
......@@ -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,
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment