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

Bézier coefficients computed for the product of two or three functions are now...

Bézier coefficients computed for the product of two or three functions are now ordered as the original convention (corner coefficients are the first)
parent f47d221b
No related branches found
No related tags found
No related merge requests found
...@@ -36,7 +36,7 @@ _CoeffDataJac::_CoeffDataJac(fullVector<double> &v, ...@@ -36,7 +36,7 @@ _CoeffDataJac::_CoeffDataJac(fullVector<double> &v,
// copy data that are not used outside of this object. // copy data that are not used outside of this object.
// It remains to swap ownership: // It remains to swap ownership:
v.setOwnData(false); v.setOwnData(false);
((fullVector<double>)_coeffs).setOwnData(true); const_cast<fullVector<double>&>(_coeffs).setOwnData(true);
_minL = _maxL = v(0); _minL = _maxL = v(0);
int i = 1; int i = 1;
...@@ -94,10 +94,8 @@ _CoeffDataScaledJac::_CoeffDataScaledJac(fullVector<double> &det, ...@@ -94,10 +94,8 @@ _CoeffDataScaledJac::_CoeffDataScaledJac(fullVector<double> &det,
// It remains to swap ownerships: // It remains to swap ownerships:
det.setOwnData(false); det.setOwnData(false);
mat.setOwnData(false); mat.setOwnData(false);
const fullVector<double> *p1 = &_coeffsJacDet; const_cast<fullVector<double>&>(_coeffsJacDet).setOwnData(true);
const fullMatrix<double> *p2 = &_coeffsJacMat; const_cast<fullMatrix<double>&>(_coeffsJacMat).setOwnData(true);
((fullVector<double>*)p1)->setOwnData(true);
((fullMatrix<double>*)p2)->setOwnData(true);
_computeAtCorner(_minL, _maxL); _computeAtCorner(_minL, _maxL);
_minB = _computeLowerBound(); _minB = _computeLowerBound();
...@@ -123,9 +121,6 @@ void _CoeffDataScaledJac::_computeAtCorner(double &min, double &max) const ...@@ -123,9 +121,6 @@ void _CoeffDataScaledJac::_computeAtCorner(double &min, double &max) const
double _CoeffDataScaledJac::_computeLowerBound() const double _CoeffDataScaledJac::_computeLowerBound() const
{ {
fullVector<double> coeffNumerator;
_bfsDet->getRaiser()->reorder(_coeffsJacDet, coeffNumerator);
fullVector<double> coeffDenominator; fullVector<double> coeffDenominator;
{ {
fullVector<double> v[3]; fullVector<double> v[3];
...@@ -145,7 +140,7 @@ double _CoeffDataScaledJac::_computeLowerBound() const ...@@ -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 bool _CoeffDataScaledJac::boundsOk(double minL, double maxL) const
...@@ -288,8 +283,6 @@ void minScaledJacobian(MElement *el, double &min) ...@@ -288,8 +283,6 @@ void minScaledJacobian(MElement *el, double &min)
min = std::min(min, domains[i]->minB()); min = std::min(min, domains[i]->minB());
delete domains[i]; delete domains[i];
} }
Msg::Info("=====================================================min %g", min);
} }
double minScaledJacobian(MElement *el) double minScaledJacobian(MElement *el)
......
...@@ -6,6 +6,22 @@ ...@@ -6,6 +6,22 @@
#include "FuncSpaceData.h" #include "FuncSpaceData.h"
#include "MElement.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) : FuncSpaceData::FuncSpaceData(const MElement *el, const bool *serendip) :
_tag(el->getTypeForMSH()), _spaceOrder(el->getPolynomialOrder()), _tag(el->getTypeForMSH()), _spaceOrder(el->getPolynomialOrder()),
_serendipity(serendip ? *serendip : el->getIsOnlySerendipity()), _serendipity(serendip ? *serendip : el->getIsOnlySerendipity()),
......
...@@ -42,6 +42,14 @@ public: ...@@ -42,6 +42,14 @@ public:
: _tag(-1), _spaceOrder(-1), _serendipity(false), _nij(-1), _nk(-1), : _tag(-1), _spaceOrder(-1), _serendipity(false), _nij(-1), _nk(-1),
_pyramidalSpace(false) {} _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* // Constructors using MElement*
FuncSpaceData(const MElement *el, const bool *serendip = NULL); FuncSpaceData(const MElement *el, const bool *serendip = NULL);
FuncSpaceData(const MElement *el, int order, const bool *serendip = NULL); FuncSpaceData(const MElement *el, int order, const bool *serendip = NULL);
...@@ -49,13 +57,12 @@ public: ...@@ -49,13 +57,12 @@ public:
bool pyr, int nij, int nk, bool pyr, int nij, int nk,
const bool *serendip = NULL); const bool *serendip = NULL);
// Constructors using element tag // Constructor using element tag
FuncSpaceData(int tag, const bool *serendip = NULL); FuncSpaceData(int tag, const bool *serendip = NULL);
// constructors using element tag or element type // constructors using element tag or element type
FuncSpaceData(bool isTag, int tagOrType, int order, FuncSpaceData(bool isTag, int tagOrType, int order,
const bool *serendip = NULL, bool elemIsSerendip = false); const bool *serendip = NULL, bool elemIsSerendip = false);
FuncSpaceData(bool isTag, int tagOrType, bool pyr, int nij, int nk, FuncSpaceData(bool isTag, int tagOrType, bool pyr, int nij, int nk,
const bool *serendip = NULL, bool elemIsSerendip = false); const bool *serendip = NULL, bool elemIsSerendip = false);
......
...@@ -442,6 +442,17 @@ fullMatrix<double> generateSubDivisorPyramid ...@@ -442,6 +442,17 @@ fullMatrix<double> generateSubDivisorPyramid
} }
return subDivisor; 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 void bezierBasis::generateBezierPoints(fullMatrix<double> &points) const
...@@ -685,25 +696,33 @@ void bezierBasisRaiser::_fillRaiserData() ...@@ -685,25 +696,33 @@ void bezierBasisRaiser::_fillRaiserData()
return; return;
} }
fullMatrix<int> exp;
{
const fullMatrix<double> &expD = _bfs->_exponents; const fullMatrix<double> &expD = _bfs->_exponents;
fullMatrix<int> exp(expD.size1(), expD.size2()); double2int(expD, exp);
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);
} }
}
int order = _bfs->getOrder(); int order = _bfs->getOrder();
int ncoeff = exp.size1(); int ncoeff = exp.size1();
int dim = _bfs->getDim(); int dim = _bfs->getDim();
int dimSimplex = _bfs->getDimSimplex(); 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; int hash = 0;
for (int l = 0; l < dim; l++) { 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++) { for (int i = 0; i < ncoeff; i++) {
...@@ -732,10 +751,29 @@ void bezierBasisRaiser::_fillRaiserData() ...@@ -732,10 +751,29 @@ void bezierBasisRaiser::_fillRaiserData()
for (int l = 0; l < dim; l++) { for (int l = 0; l < dim; l++) {
hash += (exp(i, l)+exp(j, l)) * pow_int(2*order+1, 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 i = 0; i < ncoeff; i++) {
for (int j = 0; j < ncoeff; j++) { for (int j = 0; j < ncoeff; j++) {
for (int k = 0; k < ncoeff; ++k) { for (int k = 0; k < ncoeff; ++k) {
...@@ -767,7 +805,7 @@ void bezierBasisRaiser::_fillRaiserData() ...@@ -767,7 +805,7 @@ void bezierBasisRaiser::_fillRaiserData()
for (int l = 0; l < dim; l++) { for (int l = 0; l < dim; l++) {
hash += (exp(i, l)+exp(j, l)+exp(k, l)) * pow_int(3*order+1, 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() ...@@ -785,23 +823,31 @@ void bezierBasisRaiser::_fillRaiserDataPyr()
return; return;
} }
fullMatrix<int> exp;
{
const fullMatrix<double> &expD = _bfs->_exponents; const fullMatrix<double> &expD = _bfs->_exponents;
fullMatrix<int> exp(expD.size1(), expD.size2()); double2int(expD, exp);
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);
} }
}
int ncoeff = exp.size1(); int ncoeff = exp.size1();
int order[3] = {fsdata.nij(), fsdata.nij(), fsdata.nk()}; 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; int hash = 0;
for (int l = 0; l < 3; l++) { 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++) { for (int i = 0; i < ncoeff; i++) {
...@@ -817,8 +863,28 @@ void bezierBasisRaiser::_fillRaiserDataPyr() ...@@ -817,8 +863,28 @@ void bezierBasisRaiser::_fillRaiserDataPyr()
for (int l = 0; l < 3; l++) { for (int l = 0; l < 3; l++) {
hash += (exp(i, l)+exp(j, l)) * pow_int(2*order[l]+1, 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 i = 0; i < ncoeff; i++) {
...@@ -836,7 +902,8 @@ void bezierBasisRaiser::_fillRaiserDataPyr() ...@@ -836,7 +902,8 @@ void bezierBasisRaiser::_fillRaiserDataPyr()
for (int l = 0; l < 3; l++) { for (int l = 0; l < 3; l++) {
hash += (exp(i, l)+exp(j, l)+exp(k, l)) * pow_int(3*order[l]+1, 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, ...@@ -846,15 +913,13 @@ void bezierBasisRaiser::computeCoeff(const fullVector<double> &coeffA,
const fullVector<double> &coeffB, const fullVector<double> &coeffB,
fullVector<double> &coeffSquare) fullVector<double> &coeffSquare)
{ {
coeffSquare.resize(_raiser2.size(), true); coeffSquare.resize(_Raiser2.size(), true);
std::map<int, std::vector<_Data> >::const_iterator it;
int ind = 0; for (unsigned int ind = 0; ind < _Raiser2.size(); ++ind) {
for (it = _raiser2.begin(); it != _raiser2.end(); ++it, ++ind) { for (unsigned int l = 0; l < _Raiser2[ind].size(); ++l) {
for (unsigned int l = 0; l < it->second.size(); ++l) { const int i = _Raiser2[ind][l].i;
const int i = it->second[l].i; const int j = _Raiser2[ind][l].j;
const int j = it->second[l].j; coeffSquare(ind) += _Raiser2[ind][l].val * coeffA(i) * coeffB(j);
coeffSquare(ind) += it->second[l].val * coeffA(i) * coeffB(j);
} }
} }
} }
...@@ -864,16 +929,14 @@ void bezierBasisRaiser::computeCoeff(const fullVector<double> &coeffA, ...@@ -864,16 +929,14 @@ void bezierBasisRaiser::computeCoeff(const fullVector<double> &coeffA,
const fullVector<double> &coeffC, const fullVector<double> &coeffC,
fullVector<double> &coeffCubic) fullVector<double> &coeffCubic)
{ {
coeffCubic.resize(_raiser3.size(), true); coeffCubic.resize(_Raiser3.size(), true);
std::map<int, std::vector<_Data> >::const_iterator it;
int ind = 0; for (unsigned int ind = 0; ind < _Raiser3.size(); ++ind) {
for (it = _raiser3.begin(); it != _raiser3.end(); ++it, ++ind) { for (unsigned int l = 0; l < _Raiser3[ind].size(); ++l) {
for (unsigned int l = 0; l < it->second.size(); ++l) { const int i = _Raiser3[ind][l].i;
const int i = it->second[l].i; const int j = _Raiser3[ind][l].j;
const int j = it->second[l].j; const int k = _Raiser3[ind][l].k;
const int k = it->second[l].k; coeffCubic(ind) += _Raiser3[ind][l].val * coeffA(i) * coeffB(j) * coeffC(k);
coeffCubic(ind) += it->second[l].val * coeffA(i) * coeffB(j) * coeffC(k);
} }
} }
} }
...@@ -882,16 +945,14 @@ void bezierBasisRaiser::computeCoeff(const fullMatrix<double> &coeffA, ...@@ -882,16 +945,14 @@ void bezierBasisRaiser::computeCoeff(const fullMatrix<double> &coeffA,
const fullMatrix<double> &coeffB, const fullMatrix<double> &coeffB,
fullMatrix<double> &coeffSquare) fullMatrix<double> &coeffSquare)
{ {
coeffSquare.resize(_raiser2.size(), coeffA.size2(), true); coeffSquare.resize(_Raiser2.size(), coeffA.size2(), true);
std::map<int, std::vector<_Data> >::const_iterator it;
int ind = 0; for (unsigned int ind = 0; ind < _Raiser2.size(); ++ind) {
for (it = _raiser2.begin(); it != _raiser2.end(); ++it, ++ind) { for (unsigned int l = 0; l < _Raiser2[ind].size(); ++l) {
for (unsigned int l = 0; l < it->second.size(); ++l) { const int i = _Raiser2[ind][l].i;
const int i = it->second[l].i; const int j = _Raiser2[ind][l].j;
const int j = it->second[l].j;
for (int ind2 = 0; ind2 < coeffA.size2(); ++ind2) { for (int ind2 = 0; ind2 < coeffA.size2(); ++ind2) {
coeffSquare(ind, ind2) += it->second[l].val * coeffA(i, ind2) coeffSquare(ind, ind2) += _Raiser2[ind][l].val * coeffA(i, ind2)
* coeffB(j, ind2); * coeffB(j, ind2);
} }
} }
...@@ -903,17 +964,15 @@ void bezierBasisRaiser::computeCoeff(const fullMatrix<double> &coeffA, ...@@ -903,17 +964,15 @@ void bezierBasisRaiser::computeCoeff(const fullMatrix<double> &coeffA,
const fullMatrix<double> &coeffC, const fullMatrix<double> &coeffC,
fullMatrix<double> &coeffCubic) fullMatrix<double> &coeffCubic)
{ {
coeffCubic.resize(_raiser3.size(), coeffA.size2(), true); coeffCubic.resize(_Raiser3.size(), coeffA.size2(), true);
std::map<int, std::vector<_Data> >::const_iterator it;
int ind = 0; for (unsigned int ind = 0; ind < _Raiser3.size(); ++ind) {
for (it = _raiser3.begin(); it != _raiser3.end(); ++it, ++ind) { for (unsigned int l = 0; l < _Raiser3[ind].size(); ++l) {
for (unsigned int l = 0; l < it->second.size(); ++l) { const int i = _Raiser3[ind][l].i;
const int i = it->second[l].i; const int j = _Raiser3[ind][l].j;
const int j = it->second[l].j; const int k = _Raiser3[ind][l].k;
const int k = it->second[l].k;
for (int ind2 = 0; ind2 < coeffA.size2(); ++ind2) { for (int ind2 = 0; ind2 < coeffA.size2(); ++ind2) {
coeffCubic(ind, ind2) += it->second[l].val * coeffA(i, ind2) coeffCubic(ind, ind2) += _Raiser3[ind][l].val * coeffA(i, ind2)
* coeffB(j, ind2) * coeffB(j, ind2)
* coeffC(k, ind2); * coeffC(k, ind2);
} }
...@@ -921,34 +980,3 @@ void bezierBasisRaiser::computeCoeff(const fullMatrix<double> &coeffA, ...@@ -921,34 +980,3 @@ void bezierBasisRaiser::computeCoeff(const fullMatrix<double> &coeffA,
} }
} }
void bezierBasisRaiser::reorder(const fullVector<double> &orig,
fullVector<double> &reordered)
{
reordered.resize(_raiser1.size(), false);
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) {
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);
}
}
}
}
...@@ -99,7 +99,7 @@ private : ...@@ -99,7 +99,7 @@ private :
_Data(double vv, int ii, int jj = -1, int kk = -1) : _Data(double vv, int ii, int jj = -1, int kk = -1) :
i(ii), j(jj), k(kk), val(vv) {} 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; const bezierBasis *_bfs;
public: public:
...@@ -107,8 +107,6 @@ public: ...@@ -107,8 +107,6 @@ public:
_fillRaiserData(); _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, void computeCoeff(const fullVector<double> &coeffA,
const fullVector<double> &coeffB, const fullVector<double> &coeffB,
fullVector<double> &coeffSquare); fullVector<double> &coeffSquare);
...@@ -124,13 +122,6 @@ public: ...@@ -124,13 +122,6 @@ public:
const fullMatrix<double> &coeffC, const fullMatrix<double> &coeffC,
fullMatrix<double> &coeffCubic); 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: private:
void _fillRaiserData(); void _fillRaiserData();
void _fillRaiserDataPyr(); void _fillRaiserDataPyr();
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment