diff --git a/Mesh/qualityMeasuresJacobian.cpp b/Mesh/qualityMeasuresJacobian.cpp
index 44355a0102453bf9e6ef389db8a22bcdde04eb41..6292fb8a6b3cc57fff274b95f7e57179c0d5ad55 100644
--- a/Mesh/qualityMeasuresJacobian.cpp
+++ b/Mesh/qualityMeasuresJacobian.cpp
@@ -36,7 +36,7 @@ _CoeffDataJac::_CoeffDataJac(fullVector<double> &v,
// copy data that are not used outside of this object.
// It remains to swap ownership:
v.setOwnData(false);
- ((fullVector<double>)_coeffs).setOwnData(true);
+ const_cast<fullVector<double>&>(_coeffs).setOwnData(true);
_minL = _maxL = v(0);
int i = 1;
@@ -94,10 +94,8 @@ _CoeffDataScaledJac::_CoeffDataScaledJac(fullVector<double> &det,
// It remains to swap ownerships:
det.setOwnData(false);
mat.setOwnData(false);
- const fullVector<double> *p1 = &_coeffsJacDet;
- const fullMatrix<double> *p2 = &_coeffsJacMat;
- ((fullVector<double>*)p1)->setOwnData(true);
- ((fullMatrix<double>*)p2)->setOwnData(true);
+ const_cast<fullVector<double>&>(_coeffsJacDet).setOwnData(true);
+ const_cast<fullMatrix<double>&>(_coeffsJacMat).setOwnData(true);
_computeAtCorner(_minL, _maxL);
_minB = _computeLowerBound();
@@ -123,9 +121,6 @@ void _CoeffDataScaledJac::_computeAtCorner(double &min, double &max) const
double _CoeffDataScaledJac::_computeLowerBound() const
{
- fullVector<double> coeffNumerator;
- _bfsDet->getRaiser()->reorder(_coeffsJacDet, coeffNumerator);
-
fullVector<double> coeffDenominator;
{
fullVector<double> v[3];
@@ -145,7 +140,7 @@ double _CoeffDataScaledJac::_computeLowerBound() const
}
}
- return _computeBoundRational(coeffNumerator, coeffDenominator, true);
+ return _computeBoundRational(_coeffsJacDet, coeffDenominator, true);
}
bool _CoeffDataScaledJac::boundsOk(double minL, double maxL) const
@@ -288,8 +283,6 @@ void minScaledJacobian(MElement *el, double &min)
min = std::min(min, domains[i]->minB());
delete domains[i];
}
-
- Msg::Info("=====================================================min %g", min);
}
double minScaledJacobian(MElement *el)
diff --git a/Numeric/FuncSpaceData.cpp b/Numeric/FuncSpaceData.cpp
index 8ad36a65943f5ed947c17f78c1699f80c260286c..f89be30a9287b6dac4d2b22451bf7ae05e5f0baa 100644
--- a/Numeric/FuncSpaceData.cpp
+++ b/Numeric/FuncSpaceData.cpp
@@ -6,6 +6,22 @@
#include "FuncSpaceData.h"
#include "MElement.h"
+FuncSpaceData::FuncSpaceData(const FuncSpaceData &fsd,
+ int order,
+ const bool *serendip) :
+ _tag(fsd._tag), _spaceOrder(order),
+ _serendipity(serendip ? *serendip : fsd._serendipity),
+ _nij(0), _nk(_spaceOrder), _pyramidalSpace(fsd._pyramidalSpace)
+{}
+
+FuncSpaceData::FuncSpaceData(const FuncSpaceData &fsd,
+ int nij, int nk,
+ const bool *serendip) :
+ _tag(fsd._tag), _spaceOrder(fsd._pyramidalSpace ? nij+nk : std::max(nij, nk)),
+ _serendipity(serendip ? *serendip : fsd._serendipity),
+ _nij(nij), _nk(nk), _pyramidalSpace(fsd._pyramidalSpace)
+{}
+
FuncSpaceData::FuncSpaceData(const MElement *el, const bool *serendip) :
_tag(el->getTypeForMSH()), _spaceOrder(el->getPolynomialOrder()),
_serendipity(serendip ? *serendip : el->getIsOnlySerendipity()),
diff --git a/Numeric/FuncSpaceData.h b/Numeric/FuncSpaceData.h
index 546deb052c35a8898a941ee86b445077f44fd327..5f0b889362e642d41b4aa569cffe2c36b57f0b3d 100644
--- a/Numeric/FuncSpaceData.h
+++ b/Numeric/FuncSpaceData.h
@@ -42,6 +42,14 @@ public:
: _tag(-1), _spaceOrder(-1), _serendipity(false), _nij(-1), _nk(-1),
_pyramidalSpace(false) {}
+ // Constructors of the function space of a different order
+ FuncSpaceData(const FuncSpaceData &fsd,
+ int order,
+ const bool *serendip = NULL);
+ FuncSpaceData(const FuncSpaceData &fsd,
+ int nij, int nk,
+ const bool *serendip = NULL);
+
// Constructors using MElement*
FuncSpaceData(const MElement *el, const bool *serendip = NULL);
FuncSpaceData(const MElement *el, int order, const bool *serendip = NULL);
@@ -49,13 +57,12 @@ public:
bool pyr, int nij, int nk,
const bool *serendip = NULL);
- // Constructors using element tag
+ // Constructor using element tag
FuncSpaceData(int tag, const bool *serendip = NULL);
// constructors using element tag or element type
FuncSpaceData(bool isTag, int tagOrType, int order,
const bool *serendip = NULL, bool elemIsSerendip = false);
-
FuncSpaceData(bool isTag, int tagOrType, bool pyr, int nij, int nk,
const bool *serendip = NULL, bool elemIsSerendip = false);
diff --git a/Numeric/bezierBasis.cpp b/Numeric/bezierBasis.cpp
index 13c2d61c694d724307fe0ad8f0c6d5d33f99b484..13d5f4605e41bb0d414ae23a3d49baba816cb8fe 100644
--- a/Numeric/bezierBasis.cpp
+++ b/Numeric/bezierBasis.cpp
@@ -442,6 +442,17 @@ fullMatrix<double> generateSubDivisorPyramid
}
return subDivisor;
}
+
+void double2int(const fullMatrix<double> &matDouble, fullMatrix<int> &matInt)
+{
+ matInt.resize(matDouble.size1(), matDouble.size2());
+ for (int i = 0; i < matDouble.size1(); ++i) {
+ for (int j = 0; j < matDouble.size2(); ++j) {
+ matInt(i, j) = static_cast<int>(matDouble(i, j) + .5);
+ }
+ }
+}
+
}
void bezierBasis::generateBezierPoints(fullMatrix<double> &points) const
@@ -685,25 +696,33 @@ void bezierBasisRaiser::_fillRaiserData()
return;
}
- const fullMatrix<double> &expD = _bfs->_exponents;
- fullMatrix<int> exp(expD.size1(), expD.size2());
- for (int i = 0; i < expD.size1(); ++i) {
- for (int j = 0; j < expD.size2(); ++j) {
- exp(i, j) = static_cast<int>(expD(i, j) + .5);
- }
+ fullMatrix<int> exp;
+ {
+ const fullMatrix<double> &expD = _bfs->_exponents;
+ double2int(expD, exp);
}
-
int order = _bfs->getOrder();
int ncoeff = exp.size1();
int dim = _bfs->getDim();
int dimSimplex = _bfs->getDimSimplex();
- for (int i = 0; i < ncoeff; i++) {
+ // Construction of raiser 2
+ fullMatrix<int> exp2;
+ {
+ fullMatrix<double> expD2;
+ FuncSpaceData dataRaiser2(_bfs->_data, 2*order);
+ gmshGenerateMonomials(dataRaiser2, expD2);
+ double2int(expD2, exp2);
+ _Raiser2.resize(exp2.size1());
+ }
+
+ std::map<int, int> hashToInd2;
+ for (int i = 0; i < exp2.size1(); ++i) {
int hash = 0;
for (int l = 0; l < dim; l++) {
- hash += exp(i, l) * pow_int(order+1, l);
+ hash += exp2(i, l) * pow_int(2*order+1, l);
}
- _raiser1[hash].push_back(_Data(1, i));
+ hashToInd2[hash] = i;
}
for (int i = 0; i < ncoeff; i++) {
@@ -732,10 +751,29 @@ void bezierBasisRaiser::_fillRaiserData()
for (int l = 0; l < dim; l++) {
hash += (exp(i, l)+exp(j, l)) * pow_int(2*order+1, l);
}
- _raiser2[hash].push_back(_Data(num/den, i, j));
+ _Raiser2[hashToInd2[hash]].push_back(_Data(num/den, i, j));
}
}
+ // Construction of raiser 3
+ fullMatrix<int> exp3;
+ {
+ fullMatrix<double> expD3;
+ FuncSpaceData dataRaiser3(_bfs->_data, 3*order);
+ gmshGenerateMonomials(dataRaiser3, expD3);
+ double2int(expD3, exp3);
+ _Raiser3.resize(exp3.size1());
+ }
+
+ std::map<int, int> hashToInd3;
+ for (int i = 0; i < exp3.size1(); ++i) {
+ int hash = 0;
+ for (int l = 0; l < dim; l++) {
+ hash += exp3(i, l) * pow_int(3*order+1, l);
+ }
+ hashToInd3[hash] = i;
+ }
+
for (int i = 0; i < ncoeff; i++) {
for (int j = 0; j < ncoeff; j++) {
for (int k = 0; k < ncoeff; ++k) {
@@ -767,7 +805,7 @@ void bezierBasisRaiser::_fillRaiserData()
for (int l = 0; l < dim; l++) {
hash += (exp(i, l)+exp(j, l)+exp(k, l)) * pow_int(3*order+1, l);
}
- _raiser3[hash].push_back(_Data(num/den, i, j, k));
+ _Raiser3[hashToInd3[hash]].push_back(_Data(num/den, i, j, k));
}
}
}
@@ -785,23 +823,31 @@ void bezierBasisRaiser::_fillRaiserDataPyr()
return;
}
- const fullMatrix<double> &expD = _bfs->_exponents;
- fullMatrix<int> exp(expD.size1(), expD.size2());
- for (int i = 0; i < expD.size1(); ++i) {
- for (int j = 0; j < expD.size2(); ++j) {
- exp(i, j) = static_cast<int>(expD(i, j) + .5);
- }
+ fullMatrix<int> exp;
+ {
+ const fullMatrix<double> &expD = _bfs->_exponents;
+ double2int(expD, exp);
}
-
int ncoeff = exp.size1();
int order[3] = {fsdata.nij(), fsdata.nij(), fsdata.nk()};
- for (int i = 0; i < ncoeff; i++) {
+ // Construction of raiser 2
+ fullMatrix<int> exp2;
+ {
+ fullMatrix<double> expD2;
+ FuncSpaceData dataRaiser2(_bfs->_data, 2*order[0], 2*order[2]);
+ gmshGenerateMonomials(dataRaiser2, expD2);
+ double2int(expD2, exp2);
+ _Raiser2.resize(exp2.size1());
+ }
+
+ std::map<int, int> hashToInd2;
+ for (int i = 0; i < exp2.size1(); ++i) {
int hash = 0;
for (int l = 0; l < 3; l++) {
- hash += exp(i, l) * pow_int(order[l]+1, l);
+ hash += exp2(i, l) * pow_int(2*order[l]+1, l);
}
- _raiser1[hash].push_back(_Data(1, i));
+ hashToInd2[hash] = i;
}
for (int i = 0; i < ncoeff; i++) {
@@ -817,10 +863,30 @@ void bezierBasisRaiser::_fillRaiserDataPyr()
for (int l = 0; l < 3; l++) {
hash += (exp(i, l)+exp(j, l)) * pow_int(2*order[l]+1, l);
}
- _raiser2[hash].push_back(_Data(num/den, i, j));
+ //_raiser2[hash].push_back(_Data(num/den, i, j));
+ _Raiser2[hashToInd2[hash]].push_back(_Data(num/den, i, j));
}
}
+ // Construction of raiser 3
+ fullMatrix<int> exp3;
+ {
+ fullMatrix<double> expD3;
+ FuncSpaceData dataRaiser3(_bfs->_data, 3*order[0], 3*order[2]);
+ gmshGenerateMonomials(dataRaiser3, expD3);
+ double2int(expD3, exp3);
+ _Raiser3.resize(exp3.size1());
+ }
+
+ std::map<int, int> hashToInd3;
+ for (int i = 0; i < exp3.size1(); ++i) {
+ int hash = 0;
+ for (int l = 0; l < 3; l++) {
+ hash += exp3(i, l) * pow_int(3*order[l]+1, l);
+ }
+ hashToInd3[hash] = i;
+ }
+
for (int i = 0; i < ncoeff; i++) {
for (int j = 0; j < ncoeff; j++) {
for (int k = 0; k < ncoeff; ++k) {
@@ -836,7 +902,8 @@ void bezierBasisRaiser::_fillRaiserDataPyr()
for (int l = 0; l < 3; l++) {
hash += (exp(i, l)+exp(j, l)+exp(k, l)) * pow_int(3*order[l]+1, l);
}
- _raiser3[hash].push_back(_Data(num/den, i, j, k));
+ //_raiser3[hash].push_back(_Data(num/den, i, j, k));
+ _Raiser3[hashToInd3[hash]].push_back(_Data(num/den, i, j, k));
}
}
}
@@ -846,15 +913,13 @@ void bezierBasisRaiser::computeCoeff(const fullVector<double> &coeffA,
const fullVector<double> &coeffB,
fullVector<double> &coeffSquare)
{
- coeffSquare.resize(_raiser2.size(), true);
- std::map<int, std::vector<_Data> >::const_iterator it;
-
- int ind = 0;
- for (it = _raiser2.begin(); it != _raiser2.end(); ++it, ++ind) {
- for (unsigned int l = 0; l < it->second.size(); ++l) {
- const int i = it->second[l].i;
- const int j = it->second[l].j;
- coeffSquare(ind) += it->second[l].val * coeffA(i) * coeffB(j);
+ coeffSquare.resize(_Raiser2.size(), true);
+
+ for (unsigned int ind = 0; ind < _Raiser2.size(); ++ind) {
+ for (unsigned int l = 0; l < _Raiser2[ind].size(); ++l) {
+ const int i = _Raiser2[ind][l].i;
+ const int j = _Raiser2[ind][l].j;
+ coeffSquare(ind) += _Raiser2[ind][l].val * coeffA(i) * coeffB(j);
}
}
}
@@ -864,16 +929,14 @@ void bezierBasisRaiser::computeCoeff(const fullVector<double> &coeffA,
const fullVector<double> &coeffC,
fullVector<double> &coeffCubic)
{
- coeffCubic.resize(_raiser3.size(), true);
- std::map<int, std::vector<_Data> >::const_iterator it;
-
- int ind = 0;
- for (it = _raiser3.begin(); it != _raiser3.end(); ++it, ++ind) {
- for (unsigned int l = 0; l < it->second.size(); ++l) {
- const int i = it->second[l].i;
- const int j = it->second[l].j;
- const int k = it->second[l].k;
- coeffCubic(ind) += it->second[l].val * coeffA(i) * coeffB(j) * coeffC(k);
+ coeffCubic.resize(_Raiser3.size(), true);
+
+ for (unsigned int ind = 0; ind < _Raiser3.size(); ++ind) {
+ for (unsigned int l = 0; l < _Raiser3[ind].size(); ++l) {
+ const int i = _Raiser3[ind][l].i;
+ const int j = _Raiser3[ind][l].j;
+ const int k = _Raiser3[ind][l].k;
+ coeffCubic(ind) += _Raiser3[ind][l].val * coeffA(i) * coeffB(j) * coeffC(k);
}
}
}
@@ -882,17 +945,15 @@ void bezierBasisRaiser::computeCoeff(const fullMatrix<double> &coeffA,
const fullMatrix<double> &coeffB,
fullMatrix<double> &coeffSquare)
{
- coeffSquare.resize(_raiser2.size(), coeffA.size2(), true);
- std::map<int, std::vector<_Data> >::const_iterator it;
-
- int ind = 0;
- for (it = _raiser2.begin(); it != _raiser2.end(); ++it, ++ind) {
- for (unsigned int l = 0; l < it->second.size(); ++l) {
- const int i = it->second[l].i;
- const int j = it->second[l].j;
+ coeffSquare.resize(_Raiser2.size(), coeffA.size2(), true);
+
+ for (unsigned int ind = 0; ind < _Raiser2.size(); ++ind) {
+ for (unsigned int l = 0; l < _Raiser2[ind].size(); ++l) {
+ const int i = _Raiser2[ind][l].i;
+ const int j = _Raiser2[ind][l].j;
for (int ind2 = 0; ind2 < coeffA.size2(); ++ind2) {
- coeffSquare(ind, ind2) += it->second[l].val * coeffA(i, ind2)
- * coeffB(j, ind2);
+ coeffSquare(ind, ind2) += _Raiser2[ind][l].val * coeffA(i, ind2)
+ * coeffB(j, ind2);
}
}
}
@@ -903,50 +964,17 @@ void bezierBasisRaiser::computeCoeff(const fullMatrix<double> &coeffA,
const fullMatrix<double> &coeffC,
fullMatrix<double> &coeffCubic)
{
- coeffCubic.resize(_raiser3.size(), coeffA.size2(), true);
- std::map<int, std::vector<_Data> >::const_iterator it;
-
- int ind = 0;
- for (it = _raiser3.begin(); it != _raiser3.end(); ++it, ++ind) {
- for (unsigned int l = 0; l < it->second.size(); ++l) {
- const int i = it->second[l].i;
- const int j = it->second[l].j;
- const int k = it->second[l].k;
- for (int ind2 = 0; ind2 < coeffA.size2(); ++ind2) {
- coeffCubic(ind, ind2) += it->second[l].val * coeffA(i, ind2)
- * coeffB(j, ind2)
- * coeffC(k, ind2);
- }
- }
- }
-}
-
-void bezierBasisRaiser::reorder(const fullVector<double> &orig,
- fullVector<double> &reordered)
-{
- reordered.resize(_raiser1.size(), false);
- std::map<int, std::vector<_Data> >::const_iterator it;
+ coeffCubic.resize(_Raiser3.size(), coeffA.size2(), true);
- int ind = 0;
- for (it = _raiser1.begin(); it != _raiser1.end(); ++it, ++ind) {
- for (unsigned int l = 0; l < it->second.size(); ++l) {
- reordered(ind) = orig(it->second[l].i);
- }
- }
-}
-
-void bezierBasisRaiser::reorder(const fullMatrix<double> &orig,
- fullMatrix<double> &reordered)
-{
- reordered.resize(_raiser1.size(), orig.size2());
- std::map<int, std::vector<_Data> >::const_iterator it;
-
- int ind = 0;
- for (it = _raiser1.begin(); it != _raiser1.end(); ++it, ++ind) {
- for (unsigned int l = 0; l < it->second.size(); ++l) {
- const int i = it->second[l].i;
- for (int ind2 = 0; ind2 < orig.size2(); ++ind2) {
- reordered(ind, ind2) = orig(i, ind2);
+ for (unsigned int ind = 0; ind < _Raiser3.size(); ++ind) {
+ for (unsigned int l = 0; l < _Raiser3[ind].size(); ++l) {
+ const int i = _Raiser3[ind][l].i;
+ const int j = _Raiser3[ind][l].j;
+ const int k = _Raiser3[ind][l].k;
+ for (int ind2 = 0; ind2 < coeffA.size2(); ++ind2) {
+ coeffCubic(ind, ind2) += _Raiser3[ind][l].val * coeffA(i, ind2)
+ * coeffB(j, ind2)
+ * coeffC(k, ind2);
}
}
}
diff --git a/Numeric/bezierBasis.h b/Numeric/bezierBasis.h
index 2daf0cebfbfdda874732083a14f09ee0eff8c826..8f3bdefc5d120e5e0caa8427f0e5ddbdd1c51403 100644
--- a/Numeric/bezierBasis.h
+++ b/Numeric/bezierBasis.h
@@ -99,7 +99,7 @@ private :
_Data(double vv, int ii, int jj = -1, int kk = -1) :
i(ii), j(jj), k(kk), val(vv) {}
};
- std::map<int, std::vector<_Data> > _raiser1, _raiser2, _raiser3;
+ std::vector<std::vector<_Data> > _Raiser2, _Raiser3;
const bezierBasis *_bfs;
public:
@@ -107,8 +107,6 @@ public:
_fillRaiserData();
};
- // Warning: Those method return a vector or a matrix of Bezier coefficients
- // that are not stocked in the same order than what is expected as input.
void computeCoeff(const fullVector<double> &coeffA,
const fullVector<double> &coeffB,
fullVector<double> &coeffSquare);
@@ -124,13 +122,6 @@ public:
const fullMatrix<double> &coeffC,
fullMatrix<double> &coeffCubic);
- // Method returning a vector/matrix whose coeff are in the order defined in
- // this class:
- void reorder(const fullVector<double> &orig,
- fullVector<double> &reordered);
- void reorder(const fullMatrix<double> &orig,
- fullMatrix<double> &reordered);
-
private:
void _fillRaiserData();
void _fillRaiserDataPyr();