diff --git a/Geo/MHexahedron.cpp b/Geo/MHexahedron.cpp
index ccc3796e99f565ff0b5cb8d50c8d4bd72cfb10b4..630daee66ad02052922226fc114ef08299613217 100644
--- a/Geo/MHexahedron.cpp
+++ b/Geo/MHexahedron.cpp
@@ -181,9 +181,7 @@ int MHexahedronN::getNumEdgesRep(bool curved)
const JacobianBasis* MHexahedron::getJacobianFuncSpace(int order) const
{
if (order == -1) return BasisFactory::getJacobianBasis(getTypeForMSH());
-
- int tag = ElementType::getTag(TYPE_HEX, order);
- return tag ? BasisFactory::getJacobianBasis(tag) : NULL;
+ return BasisFactory::getJacobianBasis(FuncSpaceData(this, order));
}
void _myGetFaceRep(MHexahedron *hex, int num, double *x, double *y, double *z,
diff --git a/Geo/MLine.cpp b/Geo/MLine.cpp
index ae11a4428389fb64a915b4efcce7c4efb0dfe440..73e4c7978263fad4fd41057d45800a435427e63c 100644
--- a/Geo/MLine.cpp
+++ b/Geo/MLine.cpp
@@ -11,13 +11,12 @@
#include "Context.h"
#include "qualityMeasures.h"
#include "decasteljau.h"
+#include "bezierBasis.h"
const JacobianBasis* MLine::getJacobianFuncSpace(int order) const
{
if (order == -1) return BasisFactory::getJacobianBasis(getTypeForMSH());
-
- int tag = ElementType::getTag(TYPE_LIN, order);
- return tag ? BasisFactory::getJacobianBasis(tag) : NULL;
+ return BasisFactory::getJacobianBasis(FuncSpaceData(this, order));
}
void MLine::getIntegrationPoints(int pOrder, int *npts, IntPt **pts)
diff --git a/Geo/MPrism.cpp b/Geo/MPrism.cpp
index 3ac5ae27e5dee2c7d766fc580e5fe2cbd4d6140a..12655cc28575610ca1e4f7c55907092e7c47d383 100644
--- a/Geo/MPrism.cpp
+++ b/Geo/MPrism.cpp
@@ -39,9 +39,7 @@ void MPrism::getIntegrationPoints(int pOrder, int *npts, IntPt **pts)
const JacobianBasis* MPrism::getJacobianFuncSpace(int order) const
{
if (order == -1) return BasisFactory::getJacobianBasis(getTypeForMSH());
-
- int tag = ElementType::getTag(TYPE_PRI, order);
- return tag ? BasisFactory::getJacobianBasis(tag) : NULL;
+ return BasisFactory::getJacobianBasis(FuncSpaceData(this, order));
}
double MPrism::getInnerRadius()
diff --git a/Geo/MPyramid.cpp b/Geo/MPyramid.cpp
index ea27a2a795fc51cd9fb53ce38f089a77df0886b8..b16f096029444411cacc7480c4ecb73f8e745fb5 100644
--- a/Geo/MPyramid.cpp
+++ b/Geo/MPyramid.cpp
@@ -36,9 +36,7 @@ void MPyramid::getIntegrationPoints(int pOrder, int *npts, IntPt **pts)
const JacobianBasis* MPyramid::getJacobianFuncSpace(int order) const
{
if (order == -1) return BasisFactory::getJacobianBasis(getTypeForMSH());
-
- int tag = ElementType::getTag(TYPE_PYR, order);
- return tag ? BasisFactory::getJacobianBasis(tag) : NULL;
+ return BasisFactory::getJacobianBasis(FuncSpaceData(this, order));
}
MPyramidN::~MPyramidN() {}
diff --git a/Geo/MQuadrangle.cpp b/Geo/MQuadrangle.cpp
index 468fd5306d75ab5a37f7e686ca6ed46b21894230..1e68dda69fd5b6cd2e8997944cadac1517b31e06 100644
--- a/Geo/MQuadrangle.cpp
+++ b/Geo/MQuadrangle.cpp
@@ -19,9 +19,7 @@
const JacobianBasis* MQuadrangle::getJacobianFuncSpace(int order) const
{
if (order == -1) return BasisFactory::getJacobianBasis(getTypeForMSH());
-
- int tag = ElementType::getTag(TYPE_QUA, order);
- return tag ? BasisFactory::getJacobianBasis(tag) : NULL;
+ return BasisFactory::getJacobianBasis(FuncSpaceData(this, order));
}
int MQuadrangleN::getNumEdgesRep(bool curved)
diff --git a/Geo/MTetrahedron.cpp b/Geo/MTetrahedron.cpp
index 9bda9f6b5d1c47673e416e3f08926731252969b3..14fd75c1a6f6483e210046165eafea035ba813c7 100644
--- a/Geo/MTetrahedron.cpp
+++ b/Geo/MTetrahedron.cpp
@@ -103,9 +103,7 @@ void MTetrahedron::xyz2uvw(double xyz[3], double uvw[3]) const
const JacobianBasis* MTetrahedron::getJacobianFuncSpace(int order) const
{
if (order == -1) return BasisFactory::getJacobianBasis(getTypeForMSH());
-
- int tag = ElementType::getTag(TYPE_TET, order);
- return tag ? BasisFactory::getJacobianBasis(tag) : NULL;
+ return BasisFactory::getJacobianBasis(FuncSpaceData(this, order));
}
int MTetrahedron10::getNumEdgesRep(bool curved){
diff --git a/Geo/MTriangle.cpp b/Geo/MTriangle.cpp
index 90ebc5d228c4f93f4d2a99f1e4ae2a37d66032d7..3b5638591b9b58b9acf56cf3c293b8561447ee20 100644
--- a/Geo/MTriangle.cpp
+++ b/Geo/MTriangle.cpp
@@ -117,9 +117,7 @@ void MTriangle::xyz2uvw(double xyz[3], double uvw[3]) const
const JacobianBasis* MTriangle::getJacobianFuncSpace(int order) const
{
if (order == -1) return BasisFactory::getJacobianBasis(getTypeForMSH());
-
- int tag = ElementType::getTag(TYPE_TRI, order);
- return tag ? BasisFactory::getJacobianBasis(tag) : NULL;
+ return BasisFactory::getJacobianBasis(FuncSpaceData(this, order));
}
int MTriangleN::getNumEdgesRep(bool curved) {
diff --git a/Mesh/CMakeLists.txt b/Mesh/CMakeLists.txt
index aba1b46fb2b466d6f1271679ec1fd49746ab62d3..7cf93c40fc481880ba36b735f5e78906e95d157d 100644
--- a/Mesh/CMakeLists.txt
+++ b/Mesh/CMakeLists.txt
@@ -30,7 +30,7 @@ meshGRegionRelocateVertex.cpp
BackgroundMeshManager.cpp
pointInsertionRTreeTools.cpp
pointInsertion.cpp
- qualityMeasures.cpp
+ qualityMeasures.cpp qualityMeasuresJacobian.cpp
BoundaryLayers.cpp
BDS.cpp
HighOrder.cpp
diff --git a/Mesh/qualityMeasuresJacobian.cpp b/Mesh/qualityMeasuresJacobian.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..e0978976583ba6a259835388cd3beac928f7779c
--- /dev/null
+++ b/Mesh/qualityMeasuresJacobian.cpp
@@ -0,0 +1,361 @@
+// Gmsh - Copyright (C) 1997-2015 C. Geuzaine, J.-F. Remacle
+//
+// See the LICENSE.txt file for license information. Please report all
+// bugs and problems to the public mailing list <gmsh@geuz.org>.
+
+#include "qualityMeasuresJacobian.h"
+#include "MElement.h"
+#include "BasisFactory.h"
+#include "bezierBasis.h"
+#include "JacobianBasis.h"
+
+
+namespace jacobianBasedQuality {
+
+bool _lessMinB::operator()(_CoeffData *cd1, _CoeffData *cd2) const
+{
+ return cd1->minB() > cd2->minB();
+}
+
+bool _lessMaxB::operator()(_CoeffData *cd1, _CoeffData *cd2) const
+{
+ return cd1->maxB() < cd2->maxB();
+}
+
+_CoeffDataJac::_CoeffDataJac(fullVector<double> &v,
+ const bezierBasis *bfs,
+ int depth)
+: _CoeffData(depth), _coeffs(v.getDataPtr(), v.size()), _bfs(bfs)
+{
+ if (!v.getOwnData()) {
+ Msg::Fatal("Cannot create an instance of _CoeffDataJac from a "
+ "fullVector that does not own its data.");
+ }
+ // _coeffs reuses the data of v, this avoid to allocate a new array and to
+ // copy data that are not used outside of this object.
+ // It remains to swap ownership:
+ v.setOwnData(false);
+ ((fullVector<double>)_coeffs).setOwnData(true);
+
+ _minL = _maxL = v(0);
+ int i = 1;
+ for (; i < bfs->getNumLagCoeff(); i++) {
+ _minL = std::min(_minL, v(i));
+ _maxL = std::max(_maxL, v(i));
+ }
+ _minB = _minL;
+ _maxB = _maxL;
+ for (; i < v.size(); i++) {
+ _minB = std::min(_minB, v(i));
+ _maxB = std::max(_maxB, v(i));
+ }
+}
+
+bool _CoeffDataJac::boundsOk(double minL, double maxL) const
+{
+ double tol = std::max(std::abs(minL), std::abs(maxL)) * 1e-3;
+ return (minL <= 0 || _minB > 0) &&
+ minL - _minB < tol &&
+ _maxB - maxL < tol;
+}
+
+void _CoeffDataJac::getSubCoeff(std::vector<_CoeffData*> &v) const
+{
+ v.clear();
+ v.reserve(_bfs->getNumDivision());
+ fullVector<double> subCoeff;
+ _bfs->subdivideBezCoeff(_coeffs, subCoeff);
+
+ int sz = _coeffs.size();
+ for (int i = 0; i < _bfs->getNumDivision(); i++) {
+ fullVector<double> coeff(sz);
+ coeff.copy(subCoeff, i * sz, sz, 0);
+ _CoeffDataJac *newData = new _CoeffDataJac(coeff, _bfs, _depth+1);
+ v.push_back(newData);
+ }
+}
+
+_CoeffDataScaledJac::_CoeffDataScaledJac(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 _coeffsJacMat 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 fullVector<double> *p1 = &_coeffsJacDet;
+ const fullMatrix<double> *p2 = &_coeffsJacMat;
+ ((fullVector<double>*)p1)->setOwnData(true);
+ ((fullMatrix<double>*)p2)->setOwnData(true);
+
+ _computeAtCorner(_minL, _maxL);
+ _minB = _computeLowerBound();
+ // computation of _maxB not implemented for now
+}
+
+void _CoeffDataScaledJac::_computeAtCorner(double &min, double &max) const
+{
+ min = std::numeric_limits<double>::infinity();
+ max = -min;
+
+ for (int i = 0; i < _bfsDet->getNumLagCoeff(); i++) {
+ double den = 1;
+ for (int j = 0; j < _coeffsJacMat.size2()/3; ++j) {
+ den *= std::sqrt(_coeffsJacMat(0+3*j, i) * _coeffsJacMat(0+3*j, i)+
+ _coeffsJacMat(1+3*j, i) * _coeffsJacMat(1+3*j, i)+
+ _coeffsJacMat(2+3*j, i) * _coeffsJacMat(2+3*j, i));
+ }
+ min = std::min(min, _coeffsJacDet(i) / den);
+ max = std::max(max, _coeffsJacDet(i) / den);
+ }
+}
+
+double _CoeffDataScaledJac::_computeLowerBound() const
+{
+ fullVector<double> coeffNumerator;
+ _bfsDet->getRaiser()->reorder(_coeffsJacDet, coeffNumerator);
+
+ fullVector<double> coeffDenominator;
+ {
+ fullVector<double> v[3];
+ for (int i = 0; i < _coeffsJacMat.size2()/3; ++i) {
+ v[i].resize(_coeffsJacMat.size1());
+ for (int j = 0; j < _coeffsJacMat.size1(); ++j) {
+ v[i](j) = std::sqrt(_coeffsJacMat(j, 0+3*i) * _coeffsJacMat(j, 0+3*i)+
+ _coeffsJacMat(j, 1+3*i) * _coeffsJacMat(j, 1+3*i)+
+ _coeffsJacMat(j, 2+3*i) * _coeffsJacMat(j, 2+3*i));
+ }
+ }
+ if (_coeffsJacMat.size2()/3 == 3) {
+ _bfsMat->getRaiser()->computeCoeff(v[0], v[1], v[2], coeffDenominator);
+ }
+ else {
+ _bfsMat->getRaiser()->computeCoeff(v[0], v[1], coeffDenominator);
+ }
+ }
+
+ return _computeBoundRational(coeffNumerator, coeffDenominator, true);
+}
+
+bool _CoeffDataScaledJac::boundsOk(double minL, double maxL) const
+{
+ 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);
+ v.push_back(newData);
+ }
+}
+
+template<typename Comp>
+void _subdivideDomainsMinOrMax(std::vector<_CoeffData*> &domains,
+ double &minL,
+ double &maxL)
+{
+ std::vector<_CoeffData*> subs;
+ make_heap(domains.begin(), domains.end(), Comp());
+ while (!domains[0]->boundsOk(minL, maxL)) {
+ _CoeffData *cd = domains[0];
+ pop_heap(domains.begin(), domains.end(), Comp());
+ domains.pop_back();
+ cd->getSubCoeff(subs);
+ delete cd;
+
+ for (unsigned int i = 0; i < subs.size(); i++) {
+ minL = std::min(minL, subs[i]->minL());
+ maxL = std::max(maxL, subs[i]->maxL());
+ domains.push_back(subs[i]);
+ push_heap(domains.begin(), domains.end(), Comp());
+ }
+ }
+}
+
+void _subdivideDomains(std::vector<_CoeffData*> &domains)
+{
+ if (domains.empty()) {
+ Msg::Warning("empty vector in Bezier subdivision, nothing to do");
+ return;
+ }
+ double minL = domains[0]->minL();
+ double maxL = domains[0]->maxL();
+ for (unsigned int i = 1; i < domains.size(); ++i) {
+ minL = std::min(minL, domains[i]->minL());
+ maxL = std::max(maxL, domains[i]->maxL());
+ }
+
+ _subdivideDomainsMinOrMax<_lessMinB>(domains, minL, maxL);
+ _subdivideDomainsMinOrMax<_lessMaxB>(domains, minL, maxL);
+}
+
+void minMaxJacobianDeterminant(MElement *el, double &min, double &max)
+{
+ const JacobianBasis *jfs = el->getJacobianFuncSpace();
+ if (!jfs) {
+ Msg::Error("Jacobian function space not implemented for type of element %d", el->getTypeForMSH());
+ return;
+ }
+
+ fullMatrix<double> nodesXYZ(jfs->getNumMapNodes(), 3);
+ el->getNodesCoord(nodesXYZ);
+
+ fullVector<double> coeffLag(jfs->getNumJacNodes());
+ fullVector<double> coeffBez(jfs->getNumJacNodes());
+ jfs->getSignedJacobian(nodesXYZ, coeffLag);
+ jfs->lag2Bez(coeffLag, coeffBez);
+
+ std::vector<_CoeffData*> domains;
+ domains.push_back(new _CoeffDataJac(coeffBez, jfs->getBezier(), 0));
+
+ _subdivideDomains(domains);
+
+ min = domains[0]->minB();
+ max = domains[0]->maxB();
+ delete domains[0];
+ for (unsigned int i = 1; i < domains.size(); ++i) {
+ min = std::min(min, domains[i]->minB());
+ max = std::max(max, domains[i]->maxB());
+ delete domains[i];
+ }
+}
+
+void minScaledJacobian(MElement *el, double &min)
+{
+ const int jacOrder = el->getPolynomialOrder() * el->getDim();
+ const JacobianBasis *jfs = el->getJacobianFuncSpace(jacOrder);
+ if (!jfs) {
+ Msg::Error("Jacobian function space not implemented for type of element %d", el->getTypeForMSH());
+ return;
+ }
+
+ fullMatrix<double> nodesXYZ(jfs->getNumMapNodes(), 3);
+ el->getNodesCoord(nodesXYZ);
+
+ fullVector<double> coeffDetLag(jfs->getNumJacNodes());
+ fullVector<double> coeffDetBez(jfs->getNumJacNodes());
+ jfs->getSignedJacobian(nodesXYZ, coeffDetLag);
+ jfs->lag2Bez(coeffDetLag, coeffDetBez);
+
+ const GradientBasis *gfs = BasisFactory::getGradientBasis(el->getTypeForMSH());
+
+ fullMatrix<double> coeffMatLag(gfs->getNumSamplingPoints(), 9);
+ fullMatrix<double> coeffMatBez(gfs->getNumSamplingPoints(), 9);
+ gfs->getAllGradientsFromNodes(nodesXYZ, coeffMatLag);
+ gfs->lag2Bez(coeffMatLag, coeffMatBez);
+
+ if (el->getDim() == 2) coeffMatBez.resize(coeffMatBez.size1(), 6);
+
+ std::vector<_CoeffData*> domains;
+ domains.push_back(
+ new _CoeffDataScaledJac(coeffDetBez, coeffMatBez, jfs->getBezier(),
+ gfs->getBezier(), 0)
+ );
+
+ _subdivideDomains(domains);
+
+ 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];
+ }
+
+ Msg::Info("=====================================================min %g", min);
+}
+
+double minScaledJacobian(MElement *el)
+{
+ double min;
+ minScaledJacobian(el, min);
+ return min;
+}
+
+double _computeBoundRational(const fullVector<double> &numerator,
+ const fullVector<double> &denominator,
+ bool lower, bool positiveDenom)
+{
+ //sleep(1);
+ //numerator.print("numerator");
+ //denominator.print("denominator");
+ if (numerator.size() != denominator.size()) {
+ Msg::Error("In order to compute a bound on a rational function, I need "
+ "vectors of the same size!");
+ return 0;
+ }
+
+ // upper and lower bound of the desired bound:
+ const double inf = std::numeric_limits<double>::infinity();
+ double upperBound = inf;
+ double lowerBound = -inf;
+
+ if (!positiveDenom) {
+ lower = lower ? false : true;
+ }
+
+ if (lower) {
+ // if lower is true, we seek: bound * den <= num
+ for (int i = 0; i < numerator.size(); ++i) {
+ if (denominator(i) == 0) {
+ if (numerator(i) < 0) return -inf;
+ }
+ else if (denominator(i) > 0) {
+ upperBound = std::min(upperBound, numerator(i)/denominator(i));
+ }
+ else {
+ lowerBound = std::max(lowerBound, numerator(i)/denominator(i));
+ }
+ }
+ if (lowerBound > upperBound)
+ return -inf;
+ else
+ return upperBound;
+ }
+ else {
+ // otherwise, we seek: bound * den >= num
+ for (int i = 0; i < numerator.size(); ++i) {
+ if (denominator(i) == 0) {
+ if (numerator(i) > 0) return inf;
+ }
+ else if (denominator(i) > 0) {
+ lowerBound = std::max(lowerBound, numerator(i)/denominator(i));
+ }
+ else {
+ upperBound = std::min(upperBound, numerator(i)/denominator(i));
+ }
+ }
+ if (lowerBound > upperBound)
+ return inf;
+ else
+ return lowerBound;
+ }
+}
+
+} // end namespace jacobianBasedQuality
diff --git a/Mesh/qualityMeasuresJacobian.h b/Mesh/qualityMeasuresJacobian.h
new file mode 100644
index 0000000000000000000000000000000000000000..242deb31dcd75fd26b33670983b927b2414f4dd9
--- /dev/null
+++ b/Mesh/qualityMeasuresJacobian.h
@@ -0,0 +1,98 @@
+// Gmsh - Copyright (C) 1997-2015 C. Geuzaine, J.-F. Remacle
+//
+// See the LICENSE.txt file for license information. Please report all
+// bugs and problems to the public mailing list <gmsh@geuz.org>.
+
+#ifndef _QUALITY_MEASURES_JACOBIAN_H_
+#define _QUALITY_MEASURES_JACOBIAN_H_
+
+#include <vector>
+#include "fullMatrix.h"
+
+class bezierBasis;
+class MElement;
+
+namespace jacobianBasedQuality {
+
+void minMaxJacobianDeterminant(MElement *el, double &min, double &max);
+void minScaledJacobian(MElement *el, double &min);
+double minScaledJacobian(MElement *el);
+
+class _CoeffData
+{
+protected:
+ double _minL, _maxL; //Extremum of Jac at corners
+ double _minB, _maxB; //Extremum of bezier coefficients
+ int _depth;
+
+public:
+ _CoeffData(int depth) : _minL(0), _maxL(0), _minB(0), _maxB(0),
+ _depth(depth) {}
+ virtual ~_CoeffData() {}
+
+ //inline int depth() const {return _depthSub;}
+ inline double minL() const {return _minL;}
+ inline double maxL() const {return _maxL;}
+ inline double minB() const {return _minB;}
+ inline double maxB() const {return _maxB;}
+
+ virtual bool boundsOk(double minL, double maxL) const = 0;
+ virtual void getSubCoeff(std::vector<_CoeffData*>&) const = 0;
+};
+
+struct _lessMinB {
+ bool operator()(_CoeffData*, _CoeffData*) const;
+};
+struct _lessMaxB {
+ bool operator()(_CoeffData*, _CoeffData*) const;
+};
+
+class _CoeffDataJac: public _CoeffData
+{
+private:
+ const fullVector<double> _coeffs;
+ const bezierBasis *_bfs;
+
+public:
+ _CoeffDataJac(fullVector<double> &v, const bezierBasis *bfs, int depth);
+ ~_CoeffDataJac() {}
+
+ bool boundsOk(double minL, double maxL) const;
+ void getSubCoeff(std::vector<_CoeffData*>&) const;
+};
+
+class _CoeffDataScaledJac: public _CoeffData
+{
+private:
+ const fullVector<double> _coeffsJacDet;
+ const fullMatrix<double> _coeffsJacMat;
+ const bezierBasis *_bfsDet, *_bfsMat;
+
+public:
+ _CoeffDataScaledJac(fullVector<double> &det,
+ fullMatrix<double> &mat,
+ const bezierBasis *bfsDet,
+ const bezierBasis *bfsMat,
+ int depth);
+ ~_CoeffDataScaledJac() {}
+
+ bool boundsOk(double minL, double maxL) const;
+ void getSubCoeff(std::vector<_CoeffData*>&) const;
+
+private:
+ void _computeAtCorner(double &min, double &max) const;
+ double _computeLowerBound() const;
+};
+
+//inline bool isValid(MElement *el) {
+// double min, max;
+// minMaxJacobianDeterminant(el, min, max);
+// return min > 0;
+//}
+
+double _computeBoundRational(const fullVector<double> &numerator,
+ const fullVector<double> &denominator,
+ bool lower, bool positiveDenom = true);
+}
+
+#endif
diff --git a/Numeric/JacobianBasis.cpp b/Numeric/JacobianBasis.cpp
index 1b3ef74eae642479b08719cb4bbd353acb985141..5978c1916a73cd0d557b3da25214292c504735d1 100644
--- a/Numeric/JacobianBasis.cpp
+++ b/Numeric/JacobianBasis.cpp
@@ -6,6 +6,7 @@
#include "JacobianBasis.h"
#include "pointsGenerators.h"
#include "nodalBasis.h"
+#include "bezierBasis.h"
#include "BasisFactory.h"
#include "Numeric.h"
#include <cmath>
@@ -150,11 +151,11 @@ inline void calcJDJ3D(double dxdX, double dxdY, double dxdZ,
dzdX, dzdY, dzdZ);
}
-}
+} //namespace
GradientBasis::GradientBasis(FuncSpaceData data)
- : _data(data) {
-
+ : _data(data)
+{
fullMatrix<double> samplingPoints;
gmshGeneratePoints(data, samplingPoints);
const int numSampPnts = samplingPoints.size1();
@@ -183,6 +184,11 @@ GradientBasis::GradientBasis(FuncSpaceData data)
gradShapeIdealMatY, gradShapeIdealMatZ);
}
+const bezierBasis* GradientBasis::getBezier() const
+{
+ return BasisFactory::getBezierBasis(_data);
+}
+
void GradientBasis::getGradientsFromNodes(const fullMatrix<double> &nodes,
fullMatrix<double> *dxyzdX,
fullMatrix<double> *dxyzdY,
@@ -193,6 +199,20 @@ void GradientBasis::getGradientsFromNodes(const fullMatrix<double> &nodes,
if (dxyzdZ) gradShapeMatZ.mult(nodes, *dxyzdZ);
}
+void GradientBasis::getAllGradientsFromNodes(const fullMatrix<double> &nodes,
+ fullMatrix<double> &dxyzdXYZ) const
+{
+ fullMatrix<double> prox;
+ prox.setAsProxy(dxyzdXYZ, 0, 3);
+ gradShapeMatX.mult(nodes, prox);
+
+ prox.setAsProxy(dxyzdXYZ, 3, 3);
+ gradShapeMatY.mult(nodes, prox);
+
+ prox.setAsProxy(dxyzdXYZ, 6, 3);
+ gradShapeMatZ.mult(nodes, prox);
+}
+
void GradientBasis::getIdealGradientsFromNodes(const fullMatrix<double> &nodes,
fullMatrix<double> *dxyzdX,
fullMatrix<double> *dxyzdY,
@@ -225,8 +245,14 @@ void GradientBasis::mapFromIdealElement(int type, double jac[3][3])
mapFromIdealElement(type, dxyzdX, dxyzdY, dxyzdZ);
}
+void GradientBasis::lag2Bez(const fullMatrix<double> &lag,
+ fullMatrix<double> &bez) const
+{
+ getBezier()->matrixLag2Bez.mult(lag,bez);
+}
+
JacobianBasis::JacobianBasis(FuncSpaceData data)
- : _bezier(NULL), _data(data), _dim(data.dimension())
+ : _data(data), _dim(data.dimension())
{
const int parentType = data.elementType();
const int primJacobianOrder = jacobianOrder(parentType, 1);
@@ -311,10 +337,7 @@ JacobianBasis::JacobianBasis(FuncSpaceData data)
const bezierBasis* JacobianBasis::getBezier() const
{
- if (!_bezier) {
- const_cast<JacobianBasis*>(this)->_bezier = BasisFactory::getBezierBasis(_data);
- }
- return _bezier;
+ return BasisFactory::getBezierBasis(_data);
}
// Computes (unit) normals to straight line element at barycenter (with norm of gradient as return value)
@@ -444,11 +467,11 @@ void JacobianBasis::getJacobianGeneral(int nJacNodes,
case 2 : {
fullMatrix<double> normal(1,3);
const double invScale = getPrimNormal2D(nodesXYZ, normal, idealNorm);
- if (invScale == 0) {
- for (int i = 0; i < nJacNodes; i++) jacobian(i) = 0;
- return;
- }
if (scaling) {
+ if (invScale == 0) {
+ for (int i = 0; i < nJacNodes; i++) jacobian(i) = 0;
+ return;
+ }
const double scale = 1./invScale;
normal(0,0) *= scale; normal(0,1) *= scale; normal(0,2) *= scale; // Faster to scale normal than afterwards
}
@@ -751,6 +774,18 @@ void JacobianBasis::getMetricMinAndGradients(const fullMatrix<double> &nodesXYZ,
}
}
+void JacobianBasis::lag2Bez(const fullVector<double> &lag,
+ fullVector<double> &bez) const
+{
+ getBezier()->matrixLag2Bez.mult(lag, bez);
+}
+
+void JacobianBasis::lag2Bez(const fullMatrix<double> &lag,
+ fullMatrix<double> &bez) const
+{
+ getBezier()->matrixLag2Bez.mult(lag, bez);
+}
+
// Research purpose (to be removed ?)
void JacobianBasis::interpolate(const fullVector<double> &jacobian,
const fullMatrix<double> &uvw,
diff --git a/Numeric/JacobianBasis.h b/Numeric/JacobianBasis.h
index 501f7abf39806c0e0f81bb61d4ca1e93e59225d1..d3bc5261acb909d5976759fe47a06f4a9fc4a50d 100644
--- a/Numeric/JacobianBasis.h
+++ b/Numeric/JacobianBasis.h
@@ -8,7 +8,8 @@
#include "fullMatrix.h"
#include "FuncSpaceData.h"
-#include "bezierBasis.h"
+
+class bezierBasis;
class GradientBasis {
public:
@@ -23,11 +24,14 @@ public:
int getNumSamplingPoints() const {return gradShapeMatX.size1();}
int getNumMapNodes() const {return gradShapeMatX.size2();}
+ const bezierBasis* getBezier() const;
void getGradientsFromNodes(const fullMatrix<double> &nodes,
fullMatrix<double> *dxyzdX,
fullMatrix<double> *dxyzdY,
fullMatrix<double> *dxyzdZ) const;
+ void getAllGradientsFromNodes(const fullMatrix<double> &nodes,
+ fullMatrix<double> &dxyzdXYZ) const;
void getIdealGradientsFromNodes(const fullMatrix<double> &nodes,
fullMatrix<double> *dxyzdX,
fullMatrix<double> *dxyzdY,
@@ -51,13 +55,13 @@ public:
fullVector<double> &gSVecY,
fullVector<double> &gSVecZ);
static void mapFromIdealElement(int type, double jac[3][3]);
+
+ void lag2Bez(const fullMatrix<double> &lag, fullMatrix<double> &bez) const;
};
class JacobianBasis {
private:
const GradientBasis *_gradBasis;
- const bezierBasis *_bezier;
-
const FuncSpaceData _data;
const int _dim;
@@ -65,7 +69,7 @@ private:
fullVector<double> primGradShapeBaryX, primGradShapeBaryY, primGradShapeBaryZ;
fullVector<double> primIdealGradShapeBaryX, primIdealGradShapeBaryY,
primIdealGradShapeBaryZ;
- fullMatrix<double> matrixPrimJac2Jac; // Lifts Lagrange basis of primary Jac. to Lagrange basis of Jac.
+ fullMatrix<double> matrixPrimJac2Jac; // Lifts Lagrange basis of primary Jac. to Lagrange basis of Jac.
int numJacNodes, numPrimJacNodes;
int numMapNodes, numPrimMapNodes;
@@ -177,16 +181,13 @@ public:
nodesXYZ, false, true, jacobian);
}
- inline void lag2Bez(const fullVector<double> &jac, fullVector<double> &bez) const {
- getBezier()->matrixLag2Bez.mult(jac,bez);
- }
- inline void lag2Bez(const fullMatrix<double> &jac, fullMatrix<double> &bez) const {
- getBezier()->matrixLag2Bez.mult(jac,bez);
- }
+ void lag2Bez(const fullVector<double> &lag, fullVector<double> &bez) const;
+ void lag2Bez(const fullMatrix<double> &lag, fullMatrix<double> &bez) const;
inline void primJac2Jac(const fullVector<double> &primJac, fullVector<double> &jac) const {
matrixPrimJac2Jac.mult(primJac,jac);
}
+ // Research purpose (to be removed ?)
void interpolate(const fullVector<double> &jacobian,
const fullMatrix<double> &uvw,
fullMatrix<double> &result, bool areBezier = false) const;
diff --git a/Numeric/bezierBasis.cpp b/Numeric/bezierBasis.cpp
index 34e916b8082c3171bd3b2a1c7bc44ff1e82042d8..13c2d61c694d724307fe0ad8f0c6d5d33f99b484 100644
--- a/Numeric/bezierBasis.cpp
+++ b/Numeric/bezierBasis.cpp
@@ -4,9 +4,10 @@
// bugs and problems to the public mailing list <gmsh@onelab.info>.
#include <algorithm>
+#include <vector>
+#include "bezierBasis.h"
#include "GmshDefines.h"
#include "GmshMessage.h"
-#include <vector>
#include "polynomialBasis.h"
#include "pyramidalBasis.h"
#include "pointsGenerators.h"
@@ -668,3 +669,286 @@ void bezierBasis::_constructPyr()
}
return;
}
+
+bezierBasisRaiser* bezierBasis::getRaiser() const
+{
+ if (!_raiser) {
+ const_cast<bezierBasis*>(this)->_raiser = new bezierBasisRaiser(this);
+ }
+ return _raiser;
+}
+
+void bezierBasisRaiser::_fillRaiserData()
+{
+ if (_bfs->getType() == TYPE_PYR) {
+ _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);
+ }
+ }
+
+ int order = _bfs->getOrder();
+ int ncoeff = exp.size1();
+ int dim = _bfs->getDim();
+ int dimSimplex = _bfs->getDimSimplex();
+
+ for (int i = 0; i < ncoeff; i++) {
+ int hash = 0;
+ for (int l = 0; l < dim; l++) {
+ hash += exp(i, l) * pow_int(order+1, l);
+ }
+ _raiser1[hash].push_back(_Data(1, i));
+ }
+
+ for (int i = 0; i < ncoeff; i++) {
+ for (int j = 0; j < ncoeff; j++) {
+ double num = 1, den = 1;
+ {
+ int compl1 = order;
+ int compl2 = order;
+ int compltot = 2*order;
+ for (int l = 0; l < dimSimplex; l++) {
+ num *= nChoosek(compl1, exp(i, l))
+ * nChoosek(compl2, exp(j, l));
+ den *= nChoosek(compltot, exp(i, l) + exp(j, l));
+ compl1 -= exp(i, l);
+ compl2 -= exp(j, l);
+ compltot -= exp(i, l) + exp(j, l);
+ }
+ for (int l = dimSimplex; l < dim; l++) {
+ num *= nChoosek(order, exp(i, l))
+ * nChoosek(order, exp(j, l));
+ den *= nChoosek(2*order, exp(i, l) + exp(j, l));
+ }
+ }
+
+ int hash = 0;
+ 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));
+ }
+ }
+
+ for (int i = 0; i < ncoeff; i++) {
+ for (int j = 0; j < ncoeff; j++) {
+ for (int k = 0; k < ncoeff; ++k) {
+ double num = 1, den = 1;
+ {
+ int compl1 = order;
+ int compl2 = order;
+ int compl3 = order;
+ int compltot = 3*order;
+ for (int l = 0; l < dimSimplex; l++) {
+ num *= nChoosek(compl1, exp(i, l))
+ * nChoosek(compl2, exp(j, l))
+ * nChoosek(compl3, exp(k, l));
+ den *= nChoosek(compltot, exp(i, l) + exp(j, l) + exp(k, l));
+ compl1 -= exp(i, l);
+ compl2 -= exp(j, l);
+ compl3 -= exp(k, l);
+ compltot -= exp(i, l) + exp(j, l) + exp(k, l);
+ }
+ for (int l = dimSimplex; l < dim; l++) {
+ num *= nChoosek(order, exp(i, l))
+ * nChoosek(order, exp(j, l))
+ * nChoosek(order, exp(k, l));
+ den *= nChoosek(3*order, exp(i, l) + exp(j, l) + exp(k, l));
+ }
+ }
+
+ int hash = 0;
+ 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));
+ }
+ }
+ }
+}
+
+void bezierBasisRaiser::_fillRaiserDataPyr()
+{
+ FuncSpaceData fsdata = _bfs->getFuncSpaceData();
+ if (fsdata.elementType() != TYPE_PYR) {
+ _fillRaiserData();
+ return;
+ }
+ if (fsdata.isPyramidalSpace()) {
+ Msg::Error("Bezier raiser not implemented for pyramidal space");
+ 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);
+ }
+ }
+
+ int ncoeff = exp.size1();
+ int order[3] = {fsdata.nij(), fsdata.nij(), fsdata.nk()};
+
+ for (int i = 0; i < ncoeff; i++) {
+ int hash = 0;
+ for (int l = 0; l < 3; l++) {
+ hash += exp(i, l) * pow_int(order[l]+1, l);
+ }
+ _raiser1[hash].push_back(_Data(1, i));
+ }
+
+ for (int i = 0; i < ncoeff; i++) {
+ for (int j = 0; j < ncoeff; j++) {
+ double num = 1, den = 1;
+ for (int l = 0; l < 3; l++) {
+ num *= nChoosek(order[l], exp(i, l))
+ * nChoosek(order[l], exp(j, l));
+ den *= nChoosek(2*order[l], exp(i, l) + exp(j, l));
+ }
+
+ int hash = 0;
+ 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));
+ }
+ }
+
+ for (int i = 0; i < ncoeff; i++) {
+ for (int j = 0; j < ncoeff; j++) {
+ for (int k = 0; k < ncoeff; ++k) {
+ double num = 1, den = 1;
+ for (int l = 0; l < 3; l++) {
+ num *= nChoosek(order[l], exp(i, l))
+ * nChoosek(order[l], exp(j, l))
+ * nChoosek(order[l], exp(k, l));
+ den *= nChoosek(3*order[l], exp(i, l) + exp(j, l) + exp(k, l));
+ }
+
+ int hash = 0;
+ 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));
+ }
+ }
+ }
+}
+
+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);
+ }
+ }
+}
+
+void bezierBasisRaiser::computeCoeff(const fullVector<double> &coeffA,
+ const fullVector<double> &coeffB,
+ 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);
+ }
+ }
+}
+
+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;
+ for (int ind2 = 0; ind2 < coeffA.size2(); ++ind2) {
+ coeffSquare(ind, ind2) += it->second[l].val * coeffA(i, ind2)
+ * coeffB(j, ind2);
+ }
+ }
+ }
+}
+
+void bezierBasisRaiser::computeCoeff(const fullMatrix<double> &coeffA,
+ const fullMatrix<double> &coeffB,
+ 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;
+
+ 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);
+ }
+ }
+ }
+}
+
diff --git a/Numeric/bezierBasis.h b/Numeric/bezierBasis.h
index 19ef6a5bdec77babd30fc306750e0d1c038dc330..5495f7a9433b4740021492172b4b52981da8c44d 100644
--- a/Numeric/bezierBasis.h
+++ b/Numeric/bezierBasis.h
@@ -9,19 +9,21 @@
#include <map>
#include <vector>
#include "fullMatrix.h"
-#include "ElementType.h"
#include "FuncSpaceData.h"
class MElement;
+class bezierBasisRaiser;
class bezierBasis {
private :
- // the 'numLagCoeff' first exponents are related to 'Lagrangian' values
+ // the 'numLagCoeff' first exponents are related to 'real' values
int _numLagCoeff;
int _numDivisions, _dimSimplex;
const FuncSpaceData _data;
+ bezierBasisRaiser *_raiser;
friend class MetricBasis;
+ friend class bezierBasisRaiser;
fullMatrix<double> _exponents;
public :
@@ -30,7 +32,7 @@ class bezierBasis {
fullMatrix<double> subDivisor;
// Constructors
- inline bezierBasis(FuncSpaceData data) : _data(data) {
+ inline bezierBasis(FuncSpaceData data) : _data(data), _raiser(NULL) {
if (_data.elementType() == TYPE_PYR)
_constructPyr();
else
@@ -39,11 +41,14 @@ class bezierBasis {
// get methods
inline int getDim() const {return _exponents.size2();}
+ inline int getType() const {return _data.elementType();}
inline int getOrder() const {return _data.spaceOrder();}
inline int getDimSimplex() const {return _dimSimplex;}
inline int getNumLagCoeff() const {return _numLagCoeff;}
inline int getNumDivision() const {return _numDivisions;}
inline int getNumSubNodes() const {return subDivisor.size1();}
+ inline FuncSpaceData getFuncSpaceData() const {return _data;}
+ bezierBasisRaiser* getRaiser() const;
// generate Bezier points
void generateBezierPoints(fullMatrix<double>&) const;
@@ -81,6 +86,59 @@ class bezierBasis {
void _FEpoints2BezPoints(fullMatrix<double>&) const;
};
+class bezierBasisRaiser {
+ // Let f, g, h be three function whose Bezier coefficients are given.
+ // This class allows to compute the Bezier coefficients of f*g and f*g*h.
+private :
+ class _Data;
+ std::map<int, std::vector<_Data> > _raiser1, _raiser2, _raiser3;
+ const bezierBasis *_bfs;
+
+ class _Data {
+ friend class bezierBasisRaiser;
+ private:
+ const int i, j, k;
+ const double val;
+ public:
+ _Data(double val, int i, int j = -1, int k = -1) :
+ i(i), j(j), k(k), val(val) {}
+ };
+
+public:
+ bezierBasisRaiser(const bezierBasis *bezier) : _bfs(bezier) {
+ _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);
+ void computeCoeff(const fullMatrix<double> &coeffA,
+ const fullMatrix<double> &coeffB,
+ fullMatrix<double> &coeffSquare);
+ void computeCoeff(const fullVector<double> &coeffA,
+ const fullVector<double> &coeffB,
+ const fullVector<double> &coeffC,
+ fullVector<double> &coeffCubic);
+ void computeCoeff(const fullMatrix<double> &coeffA,
+ const fullMatrix<double> &coeffB,
+ 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();
+};
+
+
#endif
diff --git a/Numeric/fullMatrix.h b/Numeric/fullMatrix.h
index bdee20667ed88fa12f12746c38f119685477c559..2bfd6f767caf849a63aa99e29a640dba919a0f19 100644
--- a/Numeric/fullMatrix.h
+++ b/Numeric/fullMatrix.h
@@ -392,6 +392,9 @@ class fullVector
{
if(fread (_data, sizeof(scalar), _r, f) != (size_t)_r) return;
}
+
+ bool getOwnData() const {return _own_data;};
+ void setOwnData(bool ownData) {_own_data = ownData;};
};
// An abstract interface for dense matrix of scalar
@@ -836,7 +839,7 @@ class fullMatrix
_data[cind+i] = x(i);
}
- bool getOwnData() {return _own_data;};
+ bool getOwnData() const {return _own_data;};
void setOwnData(bool ownData) {_own_data = ownData;};
};
#endif