diff --git a/Numeric/BasisFactory.cpp b/Numeric/BasisFactory.cpp
index 298371ee08d59cab0b2fe9c461bef1cd7ecba23d..d39b63f57cc0050abd138748d70da6c3da6ded64 100644
--- a/Numeric/BasisFactory.cpp
+++ b/Numeric/BasisFactory.cpp
@@ -14,6 +14,7 @@
std::map<int, nodalBasis*> BasisFactory::fs;
std::map<int, JacobianBasis*> BasisFactory::js;
+std::map<int, MetricBasis*> BasisFactory::ms;
BasisFactory::Cont_bezierBasis BasisFactory::bs;
BasisFactory::Cont_gradBasis BasisFactory::gs;
@@ -77,6 +78,17 @@ const JacobianBasis* BasisFactory::getJacobianBasis(int tag)
return J;
}
+const MetricBasis* BasisFactory::getMetricBasis(int tag)
+{
+ std::map<int, MetricBasis*>::const_iterator it = ms.find(tag);
+ if (it != ms.end())
+ return it->second;
+
+ MetricBasis* M = new MetricBasis(tag);
+ ms.insert(std::make_pair(tag, M));
+ return M;
+}
+
const GradientBasis* BasisFactory::getGradientBasis(int tag, int order)
{
std::pair<int, int> key(tag, order);
diff --git a/Numeric/BasisFactory.h b/Numeric/BasisFactory.h
index 171f740c340f36d1001cad3b4217561868d59ad3..6c4e4e110cf691ebe5d464a0f98e2c04482593ce 100644
--- a/Numeric/BasisFactory.h
+++ b/Numeric/BasisFactory.h
@@ -10,6 +10,7 @@
#include "MPyramid.h"
#include "nodalBasis.h"
#include "JacobianBasis.h"
+#include "MetricBasis.h"
class BasisFactory
{
@@ -19,6 +20,7 @@ class BasisFactory
private:
static std::map<int, nodalBasis*> fs;
static std::map<int, JacobianBasis*> js;
+ static std::map<int, MetricBasis*> ms;
static Cont_gradBasis gs;
static Cont_bezierBasis bs;
// store bezier bases by parentType and order (no serendipity..)
@@ -27,6 +29,8 @@ class BasisFactory
// Caution: the returned pointer can be NULL
static const nodalBasis* getNodalBasis(int tag);
static const JacobianBasis* getJacobianBasis(int tag);
+ static const MetricBasis* getMetricBasis(int tag);
+
static const GradientBasis* getGradientBasis(int tag, int order);
static const bezierBasis* getBezierBasis(int parentTag, int order);
static inline const bezierBasis* getBezierBasis(int tag) {
diff --git a/Numeric/JacobianBasis.cpp b/Numeric/JacobianBasis.cpp
index 62a4a0e778cb1d04f6ad9b125a55b02875c3a532..e23bb8ca36893c3e27357e0a2cd2001a2c7d24ff 100644
--- a/Numeric/JacobianBasis.cpp
+++ b/Numeric/JacobianBasis.cpp
@@ -659,12 +659,18 @@ void JacobianBasis::getMetricMinAndGradients(const fullMatrix<double> &nodesXYZ,
void JacobianBasis::interpolate(const fullVector<double> &jacobian,
const fullMatrix<double> &uvw,
- fullMatrix<double> &result) const
+ fullMatrix<double> &result,
+ bool areBezier) const
{
fullMatrix<double> bezM(jacobian.size(), 1);
fullVector<double> bez;
bez.setAsProxy(bezM, 0);
- lag2Bez(jacobian, bez);
+
+ if (areBezier)
+ bez.setAll(jacobian);
+ else
+ lag2Bez(jacobian, bez);
+
bezier->interpolate(bezM, uvw, result);
}
diff --git a/Numeric/JacobianBasis.h b/Numeric/JacobianBasis.h
index aec12be7e6996f1bec5fe4d9825afa0b674ab7c9..9701cec6357409eadf3974a23c10f51df599550d 100644
--- a/Numeric/JacobianBasis.h
+++ b/Numeric/JacobianBasis.h
@@ -113,7 +113,7 @@ class JacobianBasis {
//
void interpolate(const fullVector<double> &jacobian,
const fullMatrix<double> &uvw,
- fullMatrix<double> &result) const;
+ fullMatrix<double> &result, bool areBezier = false) const;
//
static int jacobianOrder(int parentType, int order);
diff --git a/Numeric/MetricBasis.cpp b/Numeric/MetricBasis.cpp
index 127b80ab69d34685db92ae88afe403f399ea608e..eadcbe3d17203dc7df2f42037a0e2d43e1c47c94 100644
--- a/Numeric/MetricBasis.cpp
+++ b/Numeric/MetricBasis.cpp
@@ -6,9 +6,14 @@
#include "MetricBasis.h"
#include "BasisFactory.h"
#include "pointsGenerators.h"
+#include "BasisFactory.h"
#include <cmath>
#include <queue>
#include "OS.h"
+#include <sstream>
+
+double MetricBasis::_tol = 1e-3;
+int MetricBasis::_which = 0;
namespace {
double cubicCardanoRoot(double p, double q)
@@ -44,11 +49,309 @@ namespace {
}
}
+MetricBasis::MetricBasis(int tag) {
+ const int type = ElementType::ParentTypeFromTag(tag);
+ const int metOrder = metricOrder(tag);
+ if (type == TYPE_HEX || type == TYPE_PRI) {
+ int order = ElementType::OrderFromTag(tag);
+ _jacobian = new JacobianBasis(tag, 3*order);
+ }
+ else if (type == TYPE_TET)
+ _jacobian = BasisFactory::getJacobianBasis(tag);
+ else
+ Msg::Fatal("metric not implemented for element tag %d", tag);
+ _gradients = BasisFactory::getGradientBasis(tag, metOrder);
+ _bezier = BasisFactory::getBezierBasis(type, metOrder);
+
+ _fillInequalities(metOrder);
+}
+
+double MetricBasis::boundRmin(const MElement *el)
+{
+ const MetricBasis *metric = BasisFactory::getMetricBasis(el->getTypeForMSH());
+ MetricData *md = NULL;
+ fullMatrix<double> dummy;
+ return metric->getBoundRmin(el, md, dummy);
+}
+
+double MetricBasis::getMinR(const MElement *el, MetricData *&md, int deg) const
+{
+ fullMatrix<double> samplingPoints;
+
+ switch (el->getType()) {
+ case TYPE_PNT :
+ samplingPoints = gmshGeneratePointsLine(0);
+ break;
+ case TYPE_LIN :
+ samplingPoints = gmshGeneratePointsLine(deg);
+ break;
+ case TYPE_TRI :
+ samplingPoints = gmshGeneratePointsTriangle(deg,false);
+ break;
+ case TYPE_QUA :
+ samplingPoints = gmshGeneratePointsQuadrangle(deg,false);
+ break;
+ case TYPE_TET :
+ samplingPoints = gmshGeneratePointsTetrahedron(deg,false);
+ break;
+ case TYPE_PRI :
+ samplingPoints = gmshGeneratePointsPrism(deg,false);
+ break;
+ case TYPE_HEX :
+ samplingPoints = gmshGeneratePointsHexahedron(deg,false);
+ break;
+ case TYPE_PYR :
+ samplingPoints = JacobianBasis::generateJacPointsPyramid(deg);
+ break;
+ default :
+ Msg::Error("Unknown Jacobian function space for element type %d", el->getType());
+ return -1;
+ }
+
+ static unsigned int aa = 0;
+ bool write = false;
+ if (md->_num < 100000 && ++aa < 200) {
+ write = true;
+ std::stringstream name;
+ name << "HoleMetric_" << el->getNum() << "_";
+ name << (md->_num % 10);
+ name << (md->_num % 100)/10;
+ name << (md->_num % 1000)/100;
+ name << (md->_num % 10000)/1000;
+ name << (md->_num % 100000)/10000;
+ name << ".txt";
+ ((MetricBasis*)this)->file.open(name.str(), std::fstream::out);
+
+ {
+ fullMatrix<double> *coeff = md->_metcoeffs;
+ fullVector<double> *jac = md->_jaccoeffs;
+ double minp, minpp, maxp, minJ2, maxJ2, minK, maxK, mina, maxa, beta, maxa2, minq, maxq, maxa3, maxK2, maxK3, RminBez, RminLag;
+ minp = _minp(*coeff);
+ minpp = _minp2(*coeff);
+ maxp = _maxp(*coeff);
+ minq = _minq(*coeff);
+ maxq = _maxq(*coeff);
+ _minMaxJacobianSqr(*jac, minJ2, maxJ2);
+ _minJ2P3(*coeff, *jac, minK);
+ _minMaxA2(*coeff, mina, maxa);
+ _maxKstA(*coeff, *jac, mina, maxK);
+
+ static const double factor2 = std::sqrt(6) * 3;
+ double x = factor2 * (minK - mina*mina*mina + mina/2);
+ double phip, term1, sin, sqrt;
+ double sq6 = std::sqrt(6);
+ {
+ if (x > 1) {
+ Msg::Warning("x = %g", x);
+ x = 1;
+ phip = M_PI / 3;
+ term1 = 1 + .5 * mina*sq6;
+ sin = std::sqrt(3) / 2;
+ sqrt = 0;
+ }
+ else if (x < -1) {
+ Msg::Warning("x = %g", x);
+ x = -1;
+ phip = 2 * M_PI / 3;
+ term1 = 1 - .5 * mina*sq6;
+ sin = std::sqrt(3) / 2;
+ sqrt = 0;
+ }
+ else {
+ phip = (std::acos(x) + M_PI) / 3;
+ term1 = 1 + mina*sq6 * std::cos(phip);
+ sin = std::sin(phip);
+ sqrt = std::sqrt(1-x*x);
+ }
+ }
+ const double dRda = sin + .5 * term1 * (1-mina*mina*sq6*sq6) / sqrt;
+ beta = -3 * mina*mina*sq6*sq6 * term1 / sqrt / dRda / 6;
+ _maxAstK3(*coeff, *jac, minK, beta, maxa2);
+ _maxAstK4(*coeff, *jac, minK, beta, maxa3);
+ _maxKstA2(*coeff, *jac, mina, beta, maxK2);
+ _maxKstA3(*coeff, *jac, mina, beta, maxK3);
+
+ if (md->_num == 22)
+ _computeRmin3(*coeff, *jac, RminLag, RminBez, 0, true);
+ else
+ _computeRmin3(*coeff, *jac, RminLag, RminBez, 0, false);
+
+ ((MetricBasis*)this)->file << minK*6*std::sqrt(6) << " ";
+ ((MetricBasis*)this)->file << maxJ2/minpp/minpp/minpp*6*std::sqrt(6) << " ";
+ ((MetricBasis*)this)->file << mina*std::sqrt(6) << " ";
+ ((MetricBasis*)this)->file << maxa*std::sqrt(6) << " ";
+ ((MetricBasis*)this)->file << beta << " " << maxa2*std::sqrt(6) << " ";
+ ((MetricBasis*)this)->file << maxK*6*std::sqrt(6) << " ";
+ ((MetricBasis*)this)->file << minp/std::sqrt(6) << " ";
+ ((MetricBasis*)this)->file << maxp/std::sqrt(6) << " ";
+ ((MetricBasis*)this)->file << minJ2 << " " << maxJ2 << " ";
+ ((MetricBasis*)this)->file << minpp/std::sqrt(6) << " ";
+ ((MetricBasis*)this)->file << minq << " " << maxq << " ";
+ ((MetricBasis*)this)->file << maxa3*std::sqrt(6) << " ";
+ ((MetricBasis*)this)->file << maxK2*6*std::sqrt(6) << " ";
+ ((MetricBasis*)this)->file << maxK3*6*std::sqrt(6) << " ";
+ ((MetricBasis*)this)->file << RminBez << " " << RminLag << std::endl;
+ }
+ }
+
+ double uvw[3];
+ double minmaxQ[2];
+ uvw[0] = samplingPoints(0, 0);
+ uvw[1] = samplingPoints(0, 1);
+ uvw[2] = samplingPoints(0, 2);
+
+ interpolate(el, md, uvw, minmaxQ, write);
+ double min, max = min = std::sqrt(minmaxQ[0]/minmaxQ[1]);
+ for (int i = 1; i < samplingPoints.size1(); ++i) {
+ uvw[0] = samplingPoints(i, 0);
+ uvw[1] = samplingPoints(i, 1);
+ uvw[2] = samplingPoints(i, 2);
+ interpolate(el, md, uvw, minmaxQ, write);
+ double tmp = std::sqrt(minmaxQ[0]/minmaxQ[1]);
+ min = std::min(min, tmp);
+ max = std::max(max, tmp);
+ //Msg::Info("%g (%g, %g)", tmp, min, max);
+ }
+ if (write) {
+ ((MetricBasis*)this)->file.close();
+ }
+ return min;
+}
+
+double MetricBasis::getBoundRmin(const MElement *el, MetricData *&md, fullMatrix<double> &lagCoeff) const
+{
+ ((MetricBasis*)this)->__curElem = (MElement*)el;
+ //if (el->getNum() != 2701) return 0;
+ int nSampPnts = _gradients->getNumSamplingPoints();
+ int nMapping = _gradients->getNumMapNodes();
+ fullMatrix<double> nodes(nMapping, 3);
+ el->getNodesCoord(nodes);
+
+ // Metric coefficients
+ fullMatrix<double> metCoeffLag;
+
+ switch (el->getDim()) {
+ case 0 :
+ return -1.;
+ case 1 :
+ case 2 :
+ Msg::Fatal("not implemented");
+ break;
+
+ case 3 :
+ {
+ fullMatrix<double> dxyzdX(nSampPnts,3), dxyzdY(nSampPnts,3), dxyzdZ(nSampPnts,3);
+ _gradients->getGradientsFromNodes(nodes, &dxyzdX, &dxyzdY, &dxyzdZ);
+
+ metCoeffLag.resize(nSampPnts, 7);
+ for (int i = 0; i < nSampPnts; i++) {
+ const double &dxdX = dxyzdX(i,0), &dydX = dxyzdX(i,1), &dzdX = dxyzdX(i,2);
+ const double &dxdY = dxyzdY(i,0), &dydY = dxyzdY(i,1), &dzdY = dxyzdY(i,2);
+ const double &dxdZ = dxyzdZ(i,0), &dydZ = dxyzdZ(i,1), &dzdZ = dxyzdZ(i,2);
+ const double dvxdX = dxdX*dxdX + dydX*dydX + dzdX*dzdX;
+ const double dvxdY = dxdY*dxdY + dydY*dydY + dzdY*dzdY;
+ const double dvxdZ = dxdZ*dxdZ + dydZ*dydZ + dzdZ*dzdZ;
+ metCoeffLag(i, 0) = (dvxdX + dvxdY + dvxdZ) / 3;
+ metCoeffLag(i, 1) = dvxdX - metCoeffLag(i, 0);
+ metCoeffLag(i, 2) = dvxdY - metCoeffLag(i, 0);
+ metCoeffLag(i, 3) = dvxdZ - metCoeffLag(i, 0);
+ const double fact = std::sqrt(2);
+ metCoeffLag(i, 4) = fact * (dxdX*dxdY + dydX*dydY + dzdX*dzdY);
+ metCoeffLag(i, 5) = fact * (dxdZ*dxdY + dydZ*dydY + dzdZ*dzdY);
+ metCoeffLag(i, 6) = fact * (dxdX*dxdZ + dydX*dydZ + dzdX*dzdZ);
+ }
+ }
+ break;
+ }
+
+ lagCoeff = metCoeffLag;
+ fullMatrix<double> *metCoeff;
+ metCoeff = new fullMatrix<double>(nSampPnts, metCoeffLag.size2());
+ _bezier->matrixLag2Bez.mult(metCoeffLag, *metCoeff);
+
+ // Jacobian coefficients
+ fullVector<double> jacLag(_jacobian->getNumJacNodes());
+ fullVector<double> *jac = new fullVector<double>(_jacobian->getNumJacNodes());
+ _jacobian->getSignedJacobian(nodes, jacLag);
+ _jacobian->lag2Bez(jacLag, *jac);
+
+ //
+ double RminLag, RminBez;
+
+ /*Msg::Info("----------------");
+ Msg::Info("Jacobian");
+ for (int i = 0; i < jac->size(); ++i) {
+ Msg::Info("%g", (*jac)(i));
+ }
+ Msg::Info("----------------");
+ Msg::Info("Metric");
+ for (int i = 0; i < metCoeff->size1(); ++i) {
+ Msg::Info("%g %g %g %g %g %g %g", (*metCoeff)(i, 0), (*metCoeff)(i, 1), (*metCoeff)(i, 2), (*metCoeff)(i, 3), (*metCoeff)(i, 4), (*metCoeff)(i, 5), (*metCoeff)(i, 6));
+ }
+ Msg::Info("----------------");*/
+
+ if (MetricBasis::_which == 1)
+ _computeRmin1(*metCoeff, *jac, RminLag, RminBez, 0);
+ else if (MetricBasis::_which == 0) {
+ _computeRmin3(*metCoeff, *jac, RminLag, RminBez, 0);
+ /*double tmpL, tmpB;
+ _computeRmin2(*metCoeff, *jac, tmpL, tmpB, 0, 3);
+ if (std::abs(tmpL-RminLag) > 1e-7) Msg::Warning("Lag %g %g", tmpL, RminLag);
+ if (std::abs(tmpB-RminBez) > 1e-3) Msg::Warning("Bez %g %g %g", tmpB, RminBez, tmpL);*/
+ }
+ else
+ _computeRmin2(*metCoeff, *jac, RminLag, RminBez, 0, MetricBasis::_which);
+
+ fullVector<double> *jjac = new fullVector<double>(*jac);
+ fullMatrix<double> *mmet = new fullMatrix<double>(*metCoeff);
+ /*for (int i = 0; i < jjac->size(); ++i) {
+ Msg::Info(":%g", (*jjac)(i));
+ }
+ for (int i = 0; i < mmet->size1(); ++i) {
+ Msg::Info(":%g | %g %g %g %g %g %g", (*mmet)(i, 0), (*mmet)(i, 1), (*mmet)(i, 2), (*mmet)(i, 3), (*mmet)(i, 4), (*mmet)(i, 5), (*mmet)(i, 6));
+ }*/
+ md = new MetricData(mmet, jjac, RminBez, 0, 0);
+ //Msg::Info("+1 %d", md);
+
+ if (RminLag-RminBez < MetricBasis::_tol) {
+ //Msg::Info("RETURNING %g", RminBez);
+ //Msg::Info("0 subdivision");
+ return RminBez;
+ }
+ else {
+ //MetricData md(metCoeff, jac, RminBez, 0);
+ MetricData *md2 = new MetricData(metCoeff, jac, RminBez, 0, 0);
+ //Msg::Info("+2 %d", md2);
+ ((MetricBasis*)this)->__numSubdivision = 0;
+ ((MetricBasis*)this)->__numSub.resize(20);
+ for (unsigned int i = 0; i < __numSub.size(); ++i) ((MetricBasis*)this)->__numSub[i] = 0;
+ ((MetricBasis*)this)->__maxdepth = 0;
+ double time = Cpu();
+ static int maxsub = 0, elmax;
+ double tt = _subdivideForRmin(md2, RminLag, MetricBasis::_tol, MetricBasis::_which);
+ if (maxsub < __numSubdivision && tt > 10-10) {
+ maxsub = __numSubdivision;
+ elmax = el->getNum();
+ }
+ Msg::Info("%d subdivisions (max %d, %d), el %d", __numSubdivision, maxsub, elmax, el->getNum());
+ /*//Msg::Info("> computation time %g", Cpu() - time);
+ Msg::Info("> maxDepth %d", __maxdepth);
+ Msg::Info("> numSubdivision %d", __numSubdivision);
+ int last = __numSub.size();
+ while (--last > 0 && __numSub[last] == 0);
+ for (unsigned int i = 0; i < last+1; ++i) {
+ Msg::Info("> depth %d: %d", i, __numSub[i]);
+ }
+ Msg::Info("RETURNING %g after subdivision", tt);*/
+ return tt;
+ }
+}
+
MetricCoefficient::MetricCoefficient(MElement *el) : _element(el)
{
const int tag = el->getTypeForMSH();
const int type = ElementType::ParentTypeFromTag(tag);
- const int metricOrder = MetricCoefficient::metricOrder(tag);
+ const int metricOrder = MetricBasis::metricOrder(tag);
if (type == TYPE_HEX || type == TYPE_PRI) {
int order = ElementType::OrderFromTag(tag);
_jacobian = new JacobianBasis(tag, 3*order);
@@ -106,11 +409,9 @@ MetricCoefficient::MetricCoefficient(MElement *el) : _element(el)
}
break;
}
-
- _fillInequalities(metricOrder);
}
-void MetricCoefficient::_fillInequalities(int metricOrder)
+void MetricBasis::_fillInequalities(int metricOrder)
{
int dimSimplex = _bezier->_dimSimplex;
int dim = _bezier->getDim();
@@ -120,10 +421,9 @@ void MetricCoefficient::_fillInequalities(int metricOrder)
exp(i, j) = static_cast<int>(_bezier->_exponents(i, j) + .5);
}
}
- int ncoeff = _coefficientsLag.size1();
+ int ncoeff = _gradients->getNumSamplingPoints();
int countP3 = 0, countJ2 = 0, countA = 0;
- double tol = .1;
for (int i = 0; i < ncoeff; i++) {
for (int j = i; j < ncoeff; j++) {
double num = 1, den = 1;
@@ -184,7 +484,9 @@ void MetricCoefficient::_fillInequalities(int metricOrder)
if (j != k) num *= 3;
}
else {
- if (j == k || i == k) num *= 3;
+ if (j == k || i == k) {
+ num *= 3;
+ }
else num *= 6;
}
@@ -193,7 +495,10 @@ void MetricCoefficient::_fillInequalities(int metricOrder)
for (int l = 0; l < dim; l++) {
hash += (exp(i, l)+exp(j, l)+exp(k, l)) * pow_int(3*metricOrder+1, l);
}
- _ineqP3[hash].push_back(IneqData(num/den, i, j, k));
+ if (j == k && j != i)
+ _ineqP3[hash].push_back(IneqData(num/den, k, j, i));
+ else
+ _ineqP3[hash].push_back(IneqData(num/den, i, j, k));
}
}
}
@@ -240,15 +545,17 @@ void MetricCoefficient::_fillInequalities(int metricOrder)
}
}
- _lightenInequalities(tol, countJ2, countP3, countA);
+ _lightenInequalities(countJ2, countP3, countA);
Msg::Info("A : %d / %d", countA, ncoeff*(ncoeff+1)/2);
Msg::Info("J2 : %d / %d", countJ2, njac*(njac+1)/2);
Msg::Info("P3 : %d / %d", countP3, ncoeff*(ncoeff+1)*(ncoeff+2)/6);
}
-void MetricCoefficient::_lightenInequalities(double tol, int &countj, int &countp, int &counta)
+
+void MetricBasis::_lightenInequalities(int &countj, int &countp, int &counta)
{
+ double tol = .0;
std::map<int, std::vector<IneqData> >::iterator it, itbeg[3], itend[3];
int cnt[3] = {0,0,0};
@@ -334,6 +641,144 @@ void MetricCoefficient::getCoefficients(fullMatrix<double> &coeff, bool bezier)
}
}
+void MetricBasis::interpolate(const MElement *el, const MetricData *md, const double *uvw, double *minmaxQ, bool write) const
+{
+ if (minmaxQ == NULL) {
+ Msg::Error("Cannot write solution of interpolation");
+ return;
+ }
+
+ int order = _bezier->getOrder();
+
+ int dimSimplex;
+ fullMatrix<double> exponents;
+ double bezuvw[3];
+ switch (el->getType()) {
+ case TYPE_PYR:
+ bezuvw[0] = .5 * (1 + uvw[0]);
+ bezuvw[1] = .5 * (1 + uvw[1]);
+ bezuvw[2] = uvw[2];
+ //_interpolateBezierPyramid(uvw, minmaxQ);
+ return;
+
+ case TYPE_HEX:
+ bezuvw[0] = .5 * (1 + uvw[0]);
+ bezuvw[1] = .5 * (1 + uvw[1]);
+ bezuvw[2] = .5 * (1 + uvw[2]);
+ dimSimplex = 0;
+ exponents = gmshGenerateMonomialsHexahedron(order);
+ break;
+
+ case TYPE_TET:
+ bezuvw[0] = uvw[0];
+ bezuvw[1] = uvw[1];
+ bezuvw[2] = uvw[2];
+ dimSimplex = 3;
+ exponents = gmshGenerateMonomialsTetrahedron(order);
+ break;
+
+ case TYPE_PRI:
+ bezuvw[0] = uvw[0];
+ bezuvw[1] = uvw[1];
+ bezuvw[2] = .5 * (1 + uvw[2]);
+ dimSimplex = 2;
+ exponents = gmshGenerateMonomialsPrism(order);
+ break;
+ }
+
+ int numCoeff = exponents.size1();
+ int dim = exponents.size2();
+
+ fullMatrix<double> metcoeffs = *md->_metcoeffs;
+ fullVector<double> jaccoeffs = *md->_jaccoeffs;
+
+ double *terms = new double[metcoeffs.size2()];
+ for (int t = 0; t < metcoeffs.size2(); ++t) {
+ terms[t] = 0;
+ for (int i = 0; i < numCoeff; i++) {
+ double dd = 1;
+ double pointCompl = 1.;
+ int exponentCompl = order;
+ for (int k = 0; k < dimSimplex; k++) {
+ dd *= nChoosek(exponentCompl, (int) exponents(i, k))
+ * pow(bezuvw[k], exponents(i, k));
+ pointCompl -= bezuvw[k];
+ exponentCompl -= (int) exponents(i, k);
+ }
+ dd *= pow(pointCompl, exponentCompl);
+
+ for (int k = dimSimplex; k < dim; k++)
+ dd *= nChoosek(order, (int) exponents(i, k))
+ * pow(bezuvw[k], exponents(i, k))
+ * pow(1. - bezuvw[k], order - exponents(i, k));
+ terms[t] += metcoeffs(i, t) * dd;
+ }
+ }
+
+ switch (metcoeffs.size2()) {
+ case 1:
+ minmaxQ[0] = terms[0];
+ minmaxQ[1] = terms[0];
+ break;
+
+ case 3:
+ {
+ double tmp = pow(terms[1], 2);
+ tmp += pow(terms[2], 2);
+ tmp = std::sqrt(tmp);
+ minmaxQ[0] = terms[0] - tmp;
+ minmaxQ[1] = terms[0] + tmp;
+ }
+ break;
+
+ case 7:
+ {
+ double tmp = pow(terms[1], 2);
+ tmp += pow(terms[2], 2);
+ tmp += pow(terms[3], 2);
+ tmp += pow(terms[4], 2);
+ tmp += pow(terms[5], 2);
+ tmp += pow(terms[6], 2);
+ tmp = std::sqrt(tmp);
+ double factor = std::sqrt(6)/3;
+ if (tmp < 1e-3*terms[0]) {
+ minmaxQ[0] = terms[0] - factor * tmp;
+ minmaxQ[1] = terms[0] + factor * tmp;
+ }
+ else {
+ double phi;
+ //{
+ fullMatrix<double> nodes(1, 3);
+ nodes(0, 0) = uvw[0];
+ nodes(0, 1) = uvw[1];
+ nodes(0, 2) = uvw[2];
+
+ fullMatrix<double> result;
+ _jacobian->interpolate(jaccoeffs, nodes, result, true);
+ phi = result(0, 0)*result(0, 0);
+ //}
+ phi -= terms[0]*terms[0]*terms[0];
+ phi += .5*terms[0]*tmp*tmp;
+ phi /= tmp*tmp*tmp;
+ phi *= 3*std::sqrt(6);
+ if (phi > 1) phi = 1;
+ if (phi < -1) phi = -1;
+ phi = std::acos(phi)/3;
+ minmaxQ[0] = terms[0] + factor * tmp * std::cos(phi + 2*M_PI/3);
+ minmaxQ[1] = terms[0] + factor * tmp * std::cos(phi);
+ ((MetricBasis*)this)->file << terms[0] << " " << tmp/std::sqrt(6) << " " << result(0, 0) << std::endl;
+ }
+ }
+ break;
+
+ default:
+ Msg::Error("Wrong number of functions for metric: %d",
+ metcoeffs.size2());
+ }
+
+ delete [] terms;
+}
+
void MetricCoefficient::interpolate(const double *uvw, double *minmaxQ)
{
if (minmaxQ == NULL) {
@@ -451,7 +896,7 @@ void MetricCoefficient::interpolate(const double *uvw, double *minmaxQ)
nodes(0, 2) = uvw[2];
fullMatrix<double> result;
- _jacobian->interpolate(jac, nodes, result);
+ _jacobian->interpolate(jac, nodes, result, false);
phi = result(0, 0)*result(0, 0);
}
phi -= terms[0]*terms[0]*terms[0];
@@ -486,9 +931,9 @@ double MetricCoefficient::getBoundRmin(double tol, int which)
double RminLag, RminBez;
if (which == 1)
- _computeRmin1(_coefficientsBez, jac, RminLag, RminBez, 0);
+ _basis->_computeRmin1(_coefficientsBez, jac, RminLag, RminBez, 0);
else
- _computeRmin2(_coefficientsBez, jac, RminLag, RminBez, 0, which);
+ _basis->_computeRmin2(_coefficientsBez, jac, RminLag, RminBez, 0, which);
if (RminLag-RminBez < tol) {
Msg::Info("RETURNING %g", RminBez);
@@ -497,13 +942,14 @@ double MetricCoefficient::getBoundRmin(double tol, int which)
else {
fullMatrix<double> *copycoeff = new fullMatrix<double>(_coefficientsBez);
fullVector<double> *copyjac = new fullVector<double>(jac);
- MetricData *md = new MetricData(copycoeff, copyjac, RminBez, 0);
+ MetricBasis::MetricData *md = new MetricBasis::MetricData(copycoeff, copyjac, RminBez, 0, 0);
+ //Msg::Info("+3 %d", md);
__numSubdivision = 0;
__numSub.resize(20);
for (unsigned int i = 0; i < __numSub.size(); ++i) __numSub[i] = 0;
__maxdepth = 0;
double time = Cpu();
- double tt = _subdivideForRmin(md, RminLag, tol, which);
+ double tt = _basis->_subdivideForRmin(md, RminLag, tol, which);
Msg::Info("> computation time %g", Cpu() - time);
Msg::Info("> maxDepth %d", __maxdepth);
Msg::Info("> numSubdivision %d", __numSubdivision);
@@ -517,7 +963,7 @@ double MetricCoefficient::getBoundRmin(double tol, int which)
}
}
-int MetricCoefficient::metricOrder(int tag)
+int MetricBasis::metricOrder(int tag)
{
const int parentType = ElementType::ParentTypeFromTag(tag);
const int order = ElementType::OrderFromTag(tag);
@@ -600,8 +1046,8 @@ void MetricCoefficient::_computeBezCoeff()
_bezier->matrixLag2Bez.mult(_coefficientsLag, _coefficientsBez);
}
-void MetricCoefficient::_computeRmin1(
- fullMatrix<double> &coeff, fullVector<double> &jac,
+void MetricBasis::_computeRmin1(
+ const fullMatrix<double> &coeff, const fullVector<double> &jac,
double &RminLag, double &RminBez, int depth) const
{
RminLag = 1.;
@@ -683,8 +1129,8 @@ void MetricCoefficient::_computeRmin1(
}
}
-void MetricCoefficient::_computeRmin2(
- fullMatrix<double> &coeff, fullVector<double> &jac,
+void MetricBasis::_computeRmin2(
+ const fullMatrix<double> &coeff, const fullVector<double> &jac,
double &RminLag, double &RminBez, int depth, int which) const
{
RminLag = 1.;
@@ -723,7 +1169,10 @@ void MetricCoefficient::_computeRmin2(
phi = std::acos(phi)/3;
p *= factor1;
double tmp = (q + p * std::cos(phi + 2*M_PI/3)) / (q + p * std::cos(phi));
- if (tmp < 0) Msg::Fatal("2 s normal ? %g", tmp);
+ if (tmp < 0) {
+ if (tmp < -1e-7) Msg::Fatal("2 s normal ? %g", tmp);
+ else tmp = 0;
+ }
RminLag = std::min(RminLag, std::sqrt(tmp));
}
}
@@ -738,9 +1187,8 @@ void MetricCoefficient::_computeRmin2(
/*if (RminBez > .6 || isnan(RminBez)) Msg::Info("minq %g maxp %g => %g", minq, maxp, RminBez);*/
}
else {
- double minJ;
+ double minJ, maxJ;
if (which == 2 || which == 4) {
- double maxJ;
_minMaxJacobianSqr(jac, minJ, maxJ);
minJ /= maxp*maxp*maxp;
}
@@ -752,7 +1200,7 @@ void MetricCoefficient::_computeRmin2(
return;
}
- double a1; //, a0;
+ double a1, a0; //TODO : a0 pas nŽcessaire
{
double p = -1./2;
double q = -minJ-1/factor2;
@@ -762,7 +1210,8 @@ void MetricCoefficient::_computeRmin2(
}
double mina;
- double maxa;
+ double maxa, maxa2;
+ double minp;
if (which == 2 || which == 5) {
mina = minq/maxp;
double minp = _minp(coeff);
@@ -773,19 +1222,14 @@ void MetricCoefficient::_computeRmin2(
}
else if (which == 3 || which == 4) {
_minMaxA2(coeff, mina, maxa);
+ _maxAstK(coeff, jac, minJ, a1, maxa2);
+ maxa = maxa2;
}
else {
Msg::Fatal("don't know what to do %d", which);
return;
}
- double phimin = std::acos(-factor1/2/mina) - M_PI/3;
- if (mina > a1) {
- RminBez = (mina+factor1*std::cos(phimin+2*M_PI/3))/(mina+factor1*std::cos(phimin)); //*
- RminBez = std::sqrt(RminBez);
- return;
- }
-
double phim = std::acos(-factor1/2/a1) - M_PI/3;
double am;
{
@@ -794,42 +1238,254 @@ void MetricCoefficient::_computeRmin2(
am = cubicCardanoRoot(p, q); //milieu
}
+
+ {
+ double pCorner[4] = {0, 0, 0, 0};
+ for (int k = 0; k < 4; ++k) {
+ for (int i = 1; i < 7; ++i) {
+ pCorner[k] += pow_int(coeff(k, i), 2);
+ }
+ pCorner[k] = std::sqrt(pCorner[k]);
+ }
+ double lagminJ, lagmina, lagmaxa;
+ lagminJ = jac(0)*jac(0) / pCorner[0]/pCorner[0]/pCorner[0];
+ lagmina = lagmaxa = coeff(0, 0) / pCorner[0];
+ static int aaa = 0;
+ //if (aaa == 0) Msg::Info(" J%g a%g", lagminJ, lagmina);
+ for (int k = 1; k < 4; ++k) {
+ lagminJ = std::min(lagminJ, jac(k)*jac(k) / pCorner[k]/pCorner[k]/pCorner[k]);
+ lagmina = std::min(lagmina, coeff(k, 0)/pCorner[k]);
+ lagmaxa = std::max(lagmaxa, coeff(k, 0)/pCorner[k]);
+ //if (aaa == 0) Msg::Info(" J%g a%g", jac(k)*jac(k) / pCorner[k]/pCorner[k]/pCorner[k], coeff(k, 0)/pCorner[k]);
+ }
+ ++aaa;
+ //Msg::Info("minJ %g[%g] (a0,am,a1)(%g, %g,%g) a(%g[%g],%g[%g],(%g)) beta%g",
+ // minJ, lagminJ, a0, am, a1, mina, lagmina, maxa, lagmaxa, maxa2, 1. / (1 - 1./6/a1/a1));
+ }
+
+
if (maxa < am) {
double phi = std::acos(factor2*(minJ-maxa*maxa*maxa+maxa/2))/3;
RminBez = (maxa+factor1*std::cos(phi+2*M_PI/3))/(maxa+factor1*std::cos(phi));
RminBez = std::sqrt(RminBez);
+ //Msg::Info("2");
return;
}
+ double phimin = std::acos(-factor1/2/mina) - M_PI/3;
if (mina > am) {
- RminBez = (mina+factor1*std::cos(phimin+2*M_PI/3))/(mina+factor1*std::cos(phimin)); //same as *
+ if (which != 2 && which != 5) {
+ minp = _minp(coeff);
+ }
+ if (which != 2 && which != 4) {
+ double tmp;
+ _minMaxJacobianSqr(jac, tmp, maxJ);
+ maxJ /= minp*minp*minp;
+ }
+ double myphi = std::acos(factor2*(maxJ-mina*mina*mina+mina/2))/3;
+ myphi = std::max(phimin, myphi);
+
+ RminBez = (mina+factor1*std::cos(myphi+2*M_PI/3))/(mina+factor1*std::cos(myphi));
+ // TODO verif plus haut si faut pas faire la mme chose du coup.
RminBez = std::sqrt(RminBez);
+ //Msg::Info("3 minJ %g maxJ %g am %g a1 %g mina %g maxa %g", minJ, maxJ, am, a1, mina, maxa);
+ //Msg::Info("3 phi %g %g %g", phimin, std::acos(factor2*(maxJ-mina*mina*mina+mina/2))/3, std::acos(-1)/3);
return;
}
else {
- RminBez = (a1+factor1*std::cos(phim+2*M_PI/3))/(a1+factor1*std::cos(phim));
+ RminBez = (am+factor1*std::cos(phim+2*M_PI/3))/(am+factor1*std::cos(phim));
RminBez = std::sqrt(RminBez);
+ //Msg::Info("4");
return;
}
}
}
-double MetricCoefficient::_subdivideForRmin(
- MetricData *md, double RminLag, double tol, int which) const
+void MetricBasis::_computeRmin3(
+ const fullMatrix<double> &coeff, const fullVector<double> &jac,
+ double &RminLag, double &RminBez,
+ int depth, bool debug) const
{
- std::priority_queue<MetricData*, std::vector<MetricData*>, lessMinB> subdomains;
- std::vector<MetricData*> vect;
- const int numCoeff = md->_metcoeffs->size2();
- const int numMetPnts = md->_metcoeffs->size1();
- const int numJacPnts = md->_jaccoeffs->size();
- const int numSub = _jacobian->getNumDivisions();
- subdomains.push(md);
+ RminLag = 1.;
+ static const double factor1 = std::sqrt(6) / 3;
+ static const double factor2 = std::sqrt(6) * 3;
- std::vector<fullVector<double>*> trash;
+ for (int i = 0; i < _bezier->getNumLagCoeff(); ++i) {
+ double q = coeff(i, 0);
+ double p = 0;
+ for (int k = 1; k < 7; ++k) {
+ p += pow_int(coeff(i, k), 2);
+ }
+ p = std::sqrt(p);
+ if (p < 1e-3 * q) {
+ p *= factor1;
+ RminLag = std::min(RminLag, std::sqrt((q - p) / (q + p))); // TODO : better
+ }
+ else {
+ double x = q/p;
+ x = .5 * x - pow_int(x,3);
+ x += jac(i)/p/p*jac(i)/p;
+ x *= factor2;
+ if (x > 1.1 || x < -1.1) {
+ if (!depth) {
+ Msg::Warning("+ phi %g (jac %g, q %g, p %g)", x, jac(i), q, p);
+ Msg::Info("%g + %g - %g = %g", jac(i)*jac(i)/p/p/p, .5 * q/p, q*q*q/p/p/p, (jac(i)*jac(i)+.5 * q*p*p-q*q*q)/p/p/p);
+ }
+ else if (depth == 1)
+ Msg::Warning("- phi %g @ %d(%d) (jac %g, q %g, p %g)", x, depth, i, jac(i), q, p);
+ }
- while (RminLag - subdomains.top()->_RminBez > tol && subdomains.size() < pow_int(8,8)) {
- fullMatrix<double> *subcoeffs, *coeff;
- fullVector<double> *subjac, *jac;
+ p *= factor1;
+ double tmpR;
+ if (x >= 1)
+ tmpR = (q - .5 * p) / (q + p);
+ else if (x <= -1)
+ tmpR = (q - p) / (q + .5 * p);
+ else {
+ const double phi = std::acos(x)/3;
+ tmpR = (q + p * std::cos(phi + 2*M_PI/3)) / (q + p * std::cos(phi));
+ }
+ //Msg::Info("%d a%g K%g R%g", i, q/p*2, jac(i)/p/p*jac(i)/p*2*2*2, tmpR);
+ if (tmpR < 0) {
+ if (tmpR < -1e-7) Msg::Fatal("3 s normal ? %g (%g, %g, %g) or (%g, %g)",
+ tmpR, p/std::sqrt(6), q, jac(i)*jac(i),
+ q/p*std::sqrt(6), jac(i)*jac(i)/p/p/p*6*std::sqrt(6));
+ else tmpR = 0;
+ }
+ //Msg::Info("tmpR %g", std::sqrt(tmpR));
+ RminLag = std::min(RminLag, std::sqrt(tmpR));
+ }
+ }
+
+ double minK;
+ _minJ2P3(coeff, jac, minK);
+ if (minK < 1e-12) {
+ RminBez = 0;
+ return;
+ }
+
+ double mina, dummy;
+ _minMaxA2(coeff, mina, dummy);
+
+ //double x = .5 * (minK - mina*mina*mina + 3*mina);
+ double x = factor2 * (minK - mina*mina*mina + mina/2);
+ double phip, term1, sin, sqrt;
+double sq6 = std::sqrt(6);
+ {
+ if (x > 1) {
+ Msg::Warning("x == %g (%g, %g)", x, mina*std::sqrt(6), minK*6*std::sqrt(6));
+ x = 1;
+ phip = M_PI / 3;
+ term1 = 1 + .5 * mina*sq6;
+ sin = std::sqrt(3) / 2;
+ sqrt = 0;
+ }
+ else if (x < -1) {
+ Msg::Warning("x == %g", x);
+ x = -1;
+ phip = 2 * M_PI / 3;
+ term1 = 1 - .5 * mina*sq6;
+ sin = std::sqrt(3) / 2;
+ sqrt = 0;
+ }
+ else {
+ phip = (std::acos(x) + M_PI) / 3;
+ term1 = 1 + mina*sq6 * std::cos(phip);
+ sin = std::sin(phip);
+ sqrt = std::sqrt(1-x*x);
+ }
+ }
+
+ const double dRda = sin + .5 * term1 * (1-mina*mina*sq6*sq6) / sqrt;
+
+ //Msg::Info("mina %g minK %g dRda %g", mina*sq6, minK*6*sq6, dRda);
+
+ //Msg::Info("a(%g, %g)", mina*sq6, dummy*sq6);
+ if (dRda < 0) {
+ //Msg::Info("dRda < 0");
+ double maxa;
+ //double maxa2;
+ double beta = -3 * mina*mina*sq6*sq6 * term1 / sqrt / dRda / 6;
+ //_maxAstK2(coeff, jac, minK, beta, maxa);
+ //_maxAstK3(coeff, jac, minK, beta, maxa2);
+ _maxAstK4(coeff, jac, minK, beta, maxa);
+ //maxa = maxa2;
+ //Msg::Info("bet%g maxa%g", beta, maxa*sq6);
+ // TODO : better am ?
+ double am, phim;
+ {
+ const double p = -1./2;
+ double q = -minK-1/factor2;
+ const double a1 = cubicCardanoRoot(p, q);
+ phim = std::acos(-factor1/2/a1) - M_PI/3;
+ q = -minK+std::cos(3*phim)/factor2;
+ am = cubicCardanoRoot(p, q);
+ }
+ if (am < maxa) {
+ RminBez = (am+factor1*std::cos(phim+2*M_PI/3))/(am+factor1*std::cos(phim));
+ if (RminBez < 0) Msg::Warning("RminBez1 = %g", RminBez);
+ RminBez = std::sqrt(RminBez);
+ return;
+ }
+ else {
+ const double phi = std::acos(factor2*(minK-maxa*maxa*maxa+maxa/2))/3;
+ RminBez = (maxa+factor1*std::cos(phi+2*M_PI/3))/(maxa+factor1*std::cos(phi));
+ if (RminBez < 0) Msg::Warning("RminBez2 = %g", RminBez);
+ RminBez = std::sqrt(RminBez);
+ return;
+ }
+ }
+ else if (term1 < 0) {
+ //double minp = _minp2(coeff);
+ //double maxJ2, dum;
+ //_minMaxJacobianSqr(jac, dum, maxJ2);
+ //double maxK = maxJ2/minp/minp/minp; // TODO enlever
+ double maxK;
+ //_maxKstA(coeff, jac, mina, maxK);
+ double beta = -3 * mina*mina*sq6*sq6 * term1 / sqrt / dRda / 6;
+ _maxKstA2(coeff, jac, mina, beta, maxK);
+ const double phimaxK = std::acos(factor2*(maxK-mina*mina*mina+mina/2))/3;
+ const double phimin = std::acos(-factor1/2/mina) - M_PI/3;
+ const double myphi = std::max(phimin, phimaxK);
+ RminBez = (mina+factor1*std::cos(myphi+2*M_PI/3))/(mina+factor1*std::cos(myphi));
+ RminBez = std::sqrt(RminBez);
+ //Msg::Info("K(%g, %g) a(%g, %g)", minK*6*sq6, maxK*6*sq6, mina*sq6, dummy*sq6);
+ if (RminBez < 0) Msg::Warning("RminBez3 = %g", RminBez);
+ return;
+ }
+ else {
+ RminBez = (mina+factor1*std::cos(phip+M_PI/3))/(mina+factor1*std::cos(phip-M_PI/3));
+ if (RminBez < 0) Msg::Warning("RminBez4 = %g", RminBez);
+ RminBez = std::sqrt(RminBez);
+ return;
+ }
+}
+
+double MetricBasis::_subdivideForRmin(
+ MetricData *md, double RminLag, double tol, int which) const
+{
+ std::priority_queue<MetricData*, std::vector<MetricData*>, lessMinB> subdomains;
+ std::vector<MetricData*> vect;
+ const int numCoeff = md->_metcoeffs->size2();
+ const int numMetPnts = md->_metcoeffs->size1();
+ const int numJacPnts = md->_jaccoeffs->size();
+ const int numSub = _jacobian->getNumDivisions();
+ subdomains.push(md);
+
+ static unsigned int aa = 0*1000;
+ bool write = false;
+ if (++aa < 200) {
+ getMinR(__curElem, md, 16);
+ }
+
+ std::vector<fullVector<double>*> trash;
+
+ //Msg::Info("lagrange %g", RminLag);
+
+ while (RminLag - subdomains.top()->_RminBez > tol && subdomains.size() < 25000) {
+ //Msg::Info("%g - %g > %g && %d < %d", RminLag, subdomains.top()->_RminBez, tol, subdomains.size(), pow_int(8,8));
+ fullMatrix<double> *subcoeffs, *coeff;
+ fullVector<double> *subjac, *jac;
MetricData *current = subdomains.top();
subcoeffs = new fullMatrix<double>(numSub*numMetPnts, numCoeff);
@@ -837,13 +1493,18 @@ double MetricCoefficient::_subdivideForRmin(
_bezier->subDivisor.mult(*current->_metcoeffs, *subcoeffs);
_jacobian->subdivideBezierCoeff(*current->_jaccoeffs, *subjac);
int depth = current->_depth;
+ int num = current->_num;
//Msg::Info("d %d RminBez %g / %g", depth, current->_RminBez, RminLag);
+
+ //Msg::Info("delete %d (%d)", current, depth);
+ //Msg::Info(" ");
delete current;
subdomains.pop();
- ++((MetricCoefficient*)this)->__numSubdivision;
- ++((MetricCoefficient*)this)->__numSub[depth];
- ((MetricCoefficient*)this)->__maxdepth = std::max(__maxdepth, depth+1);
+ ++((MetricBasis*)this)->__numSubdivision;
+ ++((MetricBasis*)this)->__numSub[depth];
+ ((MetricBasis*)this)->__maxdepth = std::max(__maxdepth, depth+1);
+ //Msg::Info("subdividing %d", current);
for (int i = 0; i < numSub; ++i) {
coeff = new fullMatrix<double>(numMetPnts, numCoeff);
@@ -853,11 +1514,21 @@ double MetricCoefficient::_subdivideForRmin(
double minLag, minBez;
if (which == 1)
_computeRmin1(*coeff, *jac, minLag, minBez, depth+1);
+ else if (which == 0)
+ _computeRmin3(*coeff, *jac, minLag, minBez, depth+1);
else
_computeRmin2(*coeff, *jac, minLag, minBez, depth+1, which);
//Msg::Info("new RminBez %g", minBez);
RminLag = std::min(RminLag, minLag);
- MetricData *metData = new MetricData(coeff, jac, minBez, depth+1);
+ int newNum = num + (i+1) * pow_int(10, depth);
+ MetricData *metData = new MetricData(coeff, jac, minBez, depth+1, newNum);
+
+ if (aa < 200) {
+ getMinR(__curElem, metData, 16);
+ }
+
+ //Msg::Info(" %g (%d)", minLag, metData);
+ //Msg::Info("+4 %d", metData);
subdomains.push(metData);
vect.push_back(metData);
}
@@ -871,23 +1542,29 @@ double MetricCoefficient::_subdivideForRmin(
return 0;*/
//Msg::Info("RminLag %g - RminBez %g @ %d", RminLag, subdomains.top()->_RminBez, subdomains.top()->_depth);
}
+ //Msg::Info("%g - %g = %g >? %g", RminLag, subdomains.top()->_RminBez, RminLag - subdomains.top()->_RminBez, tol);
+ //Msg::Info("%d <? %d", subdomains.size(), 25000);
md = subdomains.top();
double ans = md->_RminBez;
+ if (isnan(ans)) Msg::Info("ISNAN %d", subdomains.size());
while (subdomains.size()) {
md = subdomains.top();
subdomains.pop();
+ //Msg::Info("del %d", md);
+ //Msg::Info(" ");
delete md;
}
for (unsigned int i = 0; i < trash.size(); ++i) {
delete trash[i];
}
+ //Msg::Info("bez%g lag%g", ans, RminLag);
return ans;
}
-double MetricCoefficient::_minp(const fullMatrix<double> &coeff) const
+double MetricBasis::_minp(const fullMatrix<double> &coeff) const
{
fullMatrix<double> minmaxCoeff(2, 6);
for (int j = 0; j < 6; ++j) {
@@ -897,12 +1574,8 @@ double MetricCoefficient::_minp(const fullMatrix<double> &coeff) const
for (int i = 1; i < coeff.size1(); ++i) {
for (int j = 0; j < 6; ++j) {
- if (minmaxCoeff(0, j) > coeff(i, j+1)) {
- minmaxCoeff(0, j) = coeff(i, j+1);
- }
- if (minmaxCoeff(1, j) < coeff(i, j+1)) {
- minmaxCoeff(1, j) = coeff(i, j+1);
- }
+ minmaxCoeff(0, j) = std::min(coeff(i, j+1), minmaxCoeff(0, j));
+ minmaxCoeff(1, j) = std::max(coeff(i, j+1), minmaxCoeff(1, j));
}
}
@@ -917,7 +1590,29 @@ double MetricCoefficient::_minp(const fullMatrix<double> &coeff) const
return std::sqrt(ans);
}
-double MetricCoefficient::_minq(const fullMatrix<double> &coeff) const
+double MetricBasis::_minp2(const fullMatrix<double> &coeff) const
+{
+ double min = 1e10;
+ std::map<int, std::vector<IneqData> >::const_iterator it = _ineqA.begin();
+ while (it != _ineqA.end()) {
+ double val = 0;
+ for (unsigned int k = 0; k < it->second.size(); ++k) {
+ const int i = it->second[k].i;
+ const int j = it->second[k].j;
+ double tmp = 0;
+ for (int l = 1; l < 7; ++l) {
+ tmp += coeff(i, l) * coeff(j, l);
+ }
+ val += it->second[k].val * tmp;
+ }
+ min = std::min(val, min);
+ ++it;
+ }
+
+ return min > 0 ? std::sqrt(min) : 0;
+}
+
+double MetricBasis::_minq(const fullMatrix<double> &coeff) const
{
double ans = coeff(0, 0);
for (int i = 1; i < coeff.size1(); ++i) {
@@ -926,22 +1621,20 @@ double MetricCoefficient::_minq(const fullMatrix<double> &coeff) const
return ans;
}
-double MetricCoefficient::_maxp(const fullMatrix<double> &coeff) const
+double MetricBasis::_maxp(const fullMatrix<double> &coeff) const
{
double ans = 0;
for (int i = 0; i < coeff.size1(); ++i) {
- double tmp = pow_int(coeff(i, 1), 2);
- tmp += pow_int(coeff(i, 2), 2);
- tmp += pow_int(coeff(i, 3), 2);
- tmp += pow_int(coeff(i, 4), 2);
- tmp += pow_int(coeff(i, 5), 2);
- tmp += pow_int(coeff(i, 6), 2);
- if (ans < tmp) ans = tmp;
+ double tmp = 0;
+ for (int j = 1; j < 7; ++j) {
+ tmp += pow_int(coeff(i, j), 2);
+ }
+ ans = std::max(ans, tmp);
}
return std::sqrt(ans);
}
-double MetricCoefficient::_maxq(const fullMatrix<double> &coeff) const
+double MetricBasis::_maxq(const fullMatrix<double> &coeff) const
{
double ans = coeff(0, 0);
for (int i = 1; i < coeff.size1(); ++i) {
@@ -950,7 +1643,7 @@ double MetricCoefficient::_maxq(const fullMatrix<double> &coeff) const
return ans;
}
-void MetricCoefficient::_minMaxA2(
+void MetricBasis::_minMaxA2(
const fullMatrix<double> &coeff, double &min, double &max) const
{
min = 1e10;
@@ -968,21 +1661,18 @@ void MetricCoefficient::_minMaxA2(
}
den += it->second[k].val * tmp;
num += it->second[k].val * coeff(i, 0) * coeff(j, 0);
- double val = num/den;
- min = std::min(val, min);
- max = std::max(val, max);
}
+ double val = num/den;
+ min = std::min(val, min);
+ max = std::max(val, max);
++it;
}
- if (min < 0)
- Msg::Fatal("minA ne devrait pas etre negatif");
-
- min = std::sqrt(min);
+ min = min > 1/6 ? std::sqrt(min) : 1/std::sqrt(6);
max = std::sqrt(max);
}
-void MetricCoefficient::_minMaxJacobianSqr(
+void MetricBasis::_minMaxJacobianSqr(
const fullVector<double> &jac, double &min, double &max) const
{
static int a = 1;
@@ -1018,7 +1708,7 @@ void MetricCoefficient::_minMaxJacobianSqr(
}
}
-void MetricCoefficient::_minJ2P3(const fullMatrix<double> &coeff,
+void MetricBasis::_minJ2P3(const fullMatrix<double> &coeff,
const fullVector<double> &jac, double &min) const
{
fullVector<double> r(coeff.size1());
@@ -1035,6 +1725,7 @@ void MetricCoefficient::_minJ2P3(const fullMatrix<double> &coeff,
itJ = _ineqJ2.begin();
itP = _ineqP3.begin();
+ if (_ineqP3.size() != _ineqJ2.size()) Msg::Fatal("sizes P3 %d, J2 %d", _ineqP3.size(), _ineqJ2.size());
//Msg::Warning("sizes %d %d", _ineqJ2.size(), _ineqP3.size());
int count = 0;
while (itJ != _ineqJ2.end() && itP != _ineqP3.end()) {
@@ -1066,4 +1757,308 @@ void MetricCoefficient::_minJ2P3(const fullMatrix<double> &coeff,
++itP;
++count;
}
+ min = std::max(min, 0.);
+}
+
+void MetricBasis::_maxAstK(const fullMatrix<double> &coeff,
+ const fullVector<double> &jac, double minK, double a1, double &maxa) const
+{
+ fullVector<double> r(coeff.size1());
+ for (int i = 0; i < coeff.size1(); ++i) {
+ r(i) = 0;
+ for (int l = 1; l < 7; ++l) {
+ r(i) += coeff(i, l) * coeff(i, l);
+ }
+ r(i) = std::sqrt(r(i));
+ }
+
+ double beta = 1. / (1 - 1./6/a1/a1);
+ beta = 10.;
+ double min = 1e10;
+
+ std::map<int, std::vector<IneqData> >::const_iterator itJ, itP;
+ itJ = _ineqJ2.begin();
+ itP = _ineqP3.begin();
+
+ while (itJ != _ineqJ2.end() && itP != _ineqP3.end()) {
+ double num = 0, den = 0;
+ for (int l = 0; l < itJ->second.size(); ++l) {
+ const int i = itJ->second[l].i;
+ const int j = itJ->second[l].j;
+ num += itJ->second[l].val * jac(i) * jac(j);
+ }
+ num *= beta;
+ for (int l = 0; l < itP->second.size(); ++l) {
+ const int i = itP->second[l].i;
+ const int j = itP->second[l].j;
+ const int k = itP->second[l].k;
+ num -= itP->second[l].val * coeff(i, 0) * coeff(j, 0) * coeff(k, 0);
+ den += itP->second[l].val * r(i) * r(j) * r(k);
+ }
+ min = std::min(min, num/den);
+ ++itJ;
+ ++itP;
+ }
+
+ maxa = std::pow(beta*minK-min, 1/3.);
+}
+
+void MetricBasis::_maxAstK2(const fullMatrix<double> &coeff,
+ const fullVector<double> &jac, double minK, double beta, double &maxa) const
+{
+ fullVector<double> r(coeff.size1());
+ for (int i = 0; i < coeff.size1(); ++i) {
+ r(i) = 0;
+ for (int l = 1; l < 7; ++l) {
+ r(i) += coeff(i, l) * coeff(i, l);
+ }
+ r(i) = std::sqrt(r(i));
+ }
+
+ double min = 1e10;
+
+ std::map<int, std::vector<IneqData> >::const_iterator itJ, itP;
+ itJ = _ineqJ2.begin();
+ itP = _ineqP3.begin();
+
+ while (itJ != _ineqJ2.end() && itP != _ineqP3.end()) {
+ double num = 0, den = 0;
+ for (int l = 0; l < itJ->second.size(); ++l) {
+ const int i = itJ->second[l].i;
+ const int j = itJ->second[l].j;
+ num += itJ->second[l].val * jac(i) * jac(j);
+ }
+ num *= beta;
+ for (int l = 0; l < itP->second.size(); ++l) {
+ const int i = itP->second[l].i;
+ const int j = itP->second[l].j;
+ const int k = itP->second[l].k;
+ num -= itP->second[l].val * coeff(i, 0) * coeff(j, 0) * coeff(k, 0);
+ den += itP->second[l].val * r(i) * r(j) * r(k);
+ }
+ min = std::min(min, num/den);
+ ++itJ;
+ ++itP;
+ }
+
+ maxa = std::pow(beta*minK-min, 1/3.);
+}
+
+
+void MetricBasis::_maxAstK3(const fullMatrix<double> &coeff,
+ const fullVector<double> &jac, double minK, double beta, double &maxa) const
+{
+ double max = 0;
+
+ std::map<int, std::vector<IneqData> >::const_iterator itJ, itP;
+ itJ = _ineqJ2.begin();
+ itP = _ineqP3.begin();
+
+ while (itJ != _ineqJ2.end() && itP != _ineqP3.end()) {
+ double num = 0;
+ for (int l = 0; l < itJ->second.size(); ++l) {
+ const int i = itJ->second[l].i;
+ const int j = itJ->second[l].j;
+ num -= itJ->second[l].val * jac(i) * jac(j);
+ }
+ num *= beta;
+ for (int l = 0; l < itP->second.size(); ++l) {
+ const int i = itP->second[l].i;
+ const int j = itP->second[l].j;
+ const int k = itP->second[l].k;
+ num += itP->second[l].val * coeff(i, 0) * coeff(j, 0) * coeff(k, 0);
+ }
+ max = std::max(max, num);
+ ++itJ;
+ ++itP;
+ }
+ double minp = _minp2(coeff);
+ max /= minp*minp*minp;
+
+ maxa = std::pow(beta*minK+max, 1/3.);
+}
+
+void MetricBasis::_maxAstK4(const fullMatrix<double> &coeff,
+ const fullVector<double> &jac, double minK, double beta, double &maxa) const
+{
+ fullVector<double> P(coeff.size1());
+ fullMatrix<double> Q(coeff.size1(), coeff.size1());
+ for (int i = 0; i < coeff.size1(); ++i) {
+ P(i) = 0;
+ for (int l = 1; l < 7; ++l) {
+ P(i) += coeff(i, l) * coeff(i, l);
+ }
+ P(i) = std::sqrt(P(i));
+ for (int j = 0; j < coeff.size1(); ++j) {
+ Q(i, j) = 0;
+ for (int l = 1; l < 7; ++l) {
+ Q(i, j) += coeff(i, l) * coeff(j, l);
+ }
+ }
+ }
+
+ double min = 1e10;
+
+ std::map<int, std::vector<IneqData> >::const_iterator itJ, itP;
+ itJ = _ineqJ2.begin();
+ itP = _ineqP3.begin();
+
+ while (itJ != _ineqJ2.end() && itP != _ineqP3.end()) {
+ double num = 0, den = 0;
+ for (int l = 0; l < itJ->second.size(); ++l) {
+ const int i = itJ->second[l].i;
+ const int j = itJ->second[l].j;
+ num += itJ->second[l].val * jac(i) * jac(j);
+ }
+ num *= beta;
+ for (int l = 0; l < itP->second.size(); ++l) {
+ const int i = itP->second[l].i;
+ const int j = itP->second[l].j;
+ const int k = itP->second[l].k;
+ num -= itP->second[l].val * coeff(i, 0) * coeff(j, 0) * coeff(k, 0);
+ double tmp = P(i) * Q(j, k);
+ tmp = std::min(min, P(j) * Q(i, k));
+ tmp = std::min(min, P(k) * Q(i, j));
+ tmp = P(i) * P(j) * P(k);
+ den += itP->second[l].val * tmp;
+ }
+ min = std::min(min, num/den);
+ ++itJ;
+ ++itP;
+ }
+
+ maxa = std::pow(beta*minK-min, 1/3.);
+}
+
+
+void MetricBasis::_maxKstA(const fullMatrix<double> &coeff,
+ const fullVector<double> &jac, double mina, double &maxK) const
+{
+ fullVector<double> r(coeff.size1());
+ for (int i = 0; i < coeff.size1(); ++i) {
+ r(i) = 0;
+ for (int l = 1; l < 7; ++l) {
+ r(i) += coeff(i, l) * coeff(i, l);
+ }
+ r(i) = std::sqrt(r(i));
+ }
+
+ double min = 1e10;
+
+ std::map<int, std::vector<IneqData> >::const_iterator itJ, itP;
+ itJ = _ineqJ2.begin();
+ itP = _ineqP3.begin();
+
+ while (itJ != _ineqJ2.end() && itP != _ineqP3.end()) {
+ double num = 0, den = 0;
+ for (int l = 0; l < itJ->second.size(); ++l) {
+ const int i = itJ->second[l].i;
+ const int j = itJ->second[l].j;
+ num -= itJ->second[l].val * jac(i) * jac(j);
+ }
+ for (int l = 0; l < itP->second.size(); ++l) {
+ const int i = itP->second[l].i;
+ const int j = itP->second[l].j;
+ const int k = itP->second[l].k;
+ num += itP->second[l].val * coeff(i, 0) * coeff(j, 0) * coeff(k, 0);
+ den += itP->second[l].val * r(i) * r(j) * r(k);
+ }
+ min = std::min(min, num/den);
+ ++itJ;
+ ++itP;
+ }
+
+ maxK = mina*mina*mina-min;
+}
+
+void MetricBasis::_maxKstA2(const fullMatrix<double> &coeff,
+ const fullVector<double> &jac, double mina, double beta, double &maxK) const
+{
+ fullVector<double> r(coeff.size1());
+ for (int i = 0; i < coeff.size1(); ++i) {
+ r(i) = 0;
+ for (int l = 1; l < 7; ++l) {
+ r(i) += coeff(i, l) * coeff(i, l);
+ }
+ r(i) = std::sqrt(r(i));
+ }
+
+ double min = 1e10;
+
+ std::map<int, std::vector<IneqData> >::const_iterator itJ, itP;
+ itJ = _ineqJ2.begin();
+ itP = _ineqP3.begin();
+
+ while (itJ != _ineqJ2.end() && itP != _ineqP3.end()) {
+ double num = 0, den = 0;
+ for (int l = 0; l < itJ->second.size(); ++l) {
+ const int i = itJ->second[l].i;
+ const int j = itJ->second[l].j;
+ num -= itJ->second[l].val * jac(i) * jac(j);
+ }
+ num *= beta;
+ for (int l = 0; l < itP->second.size(); ++l) {
+ const int i = itP->second[l].i;
+ const int j = itP->second[l].j;
+ const int k = itP->second[l].k;
+ num += itP->second[l].val * coeff(i, 0) * coeff(j, 0) * coeff(k, 0);
+ den += itP->second[l].val * r(i) * r(j) * r(k);
+ }
+ min = std::min(min, num/den);
+ ++itJ;
+ ++itP;
+ }
+
+ maxK = 1/beta*(mina*mina*mina-min);
+}
+
+void MetricBasis::_maxKstA3(const fullMatrix<double> &coeff,
+ const fullVector<double> &jac, double mina, double beta, double &maxK) const
+{
+ fullVector<double> P(coeff.size1());
+ fullMatrix<double> Q(coeff.size1(), coeff.size1());
+ for (int i = 0; i < coeff.size1(); ++i) {
+ P(i) = 0;
+ for (int l = 1; l < 7; ++l) {
+ P(i) += coeff(i, l) * coeff(i, l);
+ }
+ P(i) = std::sqrt(P(i));
+ for (int j = 0; j < coeff.size1(); ++j) {
+ Q(i, j) = 0;
+ for (int l = 1; l < 7; ++l) {
+ Q(i, j) += coeff(i, l) * coeff(j, l);
+ }
+ }
+ }
+
+ double min = 1e10;
+
+ std::map<int, std::vector<IneqData> >::const_iterator itJ, itP;
+ itJ = _ineqJ2.begin();
+ itP = _ineqP3.begin();
+
+ while (itJ != _ineqJ2.end() && itP != _ineqP3.end()) {
+ double num = 0, den = 0;
+ for (int l = 0; l < itJ->second.size(); ++l) {
+ const int i = itJ->second[l].i;
+ const int j = itJ->second[l].j;
+ num -= itJ->second[l].val * jac(i) * jac(j);
+ }
+ num *= beta;
+ for (int l = 0; l < itP->second.size(); ++l) {
+ const int i = itP->second[l].i;
+ const int j = itP->second[l].j;
+ const int k = itP->second[l].k;
+ num += itP->second[l].val * coeff(i, 0) * coeff(j, 0) * coeff(k, 0);
+ if (j == k)
+ den += itP->second[l].val * Q(i,i) * P(i);
+ else
+ den += itP->second[l].val * 1/3*(Q(i,j)*P(k)+Q(i,k)*P(j)+Q(k,j)*P(i));
+ }
+ min = std::min(min, num/den);
+ ++itJ;
+ ++itP;
+ }
+
+ maxK = 1/beta*(mina*mina*mina-min);
}
diff --git a/Numeric/MetricBasis.h b/Numeric/MetricBasis.h
index eadbfd8f49aade6e489d94a7af8cafc969f4328d..b57effcf887219f12150e72e531766b0cc62af90 100644
--- a/Numeric/MetricBasis.h
+++ b/Numeric/MetricBasis.h
@@ -9,42 +9,138 @@
#include "MElement.h"
#include "JacobianBasis.h"
#include "fullMatrix.h"
+#include <fstream>
-class MetricCoefficient {
+class MetricBasis {
+ friend class MetricCoefficient;
+ friend class GMSH_AnalyseCurvedMeshPlugin;
private:
+ const JacobianBasis *_jacobian;
+ const GradientBasis *_gradients;
+ const bezierBasis *_bezier;
+ static double _tol;
+ static int _which;
+
+ int __maxdepth, __numSubdivision;
+ std::vector<int> __numSub;
+ MElement *__curElem;
+
+ std::fstream file;
+
+ class IneqData {
+ public:
+ int i, j, k;
+ double val;
+
+ public:
+ IneqData(double val, int i, int j, int k = -1) : i(i), j(j), k(k), val(val) {}
+ };
+
class MetricData {
public:
fullMatrix<double> *_metcoeffs;
fullVector<double> *_jaccoeffs;
double _RminBez;
- int _depth;
+ int _depth, _num;
public:
- MetricData(fullMatrix<double> *m, fullVector<double> *j, double r, int d) :
- _metcoeffs(m), _jaccoeffs(j), _RminBez(r), _depth(d) {}
+ MetricData(fullMatrix<double> *m, fullVector<double> *j, double r, int d, int num) :
+ _metcoeffs(m), _jaccoeffs(j), _RminBez(r), _depth(d), _num(num) {}
~MetricData() {
delete _metcoeffs;
delete _jaccoeffs;
}
};
+
+ std::map<int, std::vector<IneqData> > _ineqJ2, _ineqP3, _ineqA;
+
+public:
+ MetricBasis(int elementTag);
+
+ static void setTol(double tol) {_tol = tol;}
+ static double getTol() {return _tol;}
+ static void setWhich(int which) {_which = which;}
+
+ double getBoundRmin(const MElement*, MetricData*&, fullMatrix<double>&) const;
+ double getMinR(const MElement*, MetricData*&, int) const;
+ static double boundRmin(const MElement *el);
+ //double getBoundRmin(int, MElement**, double*);
+ //static double boundRmin(int, MElement**, double*, bool sameType = false);
+
+ void interpolate(const MElement*, const MetricData*, const double *uvw, double *minmaxQ, bool write = false) const;
+
+ static int metricOrder(int tag);
+
+private:
+ void _fillInequalities(int order);
+ void _lightenInequalities(int&, int&, int&); //TODO change
+
+ void _computeRmin1(const fullMatrix<double>&, const fullVector<double>&,
+ double &RminLag, double &RminBez, int) const;
+ void _computeRmin2(const fullMatrix<double>&, const fullVector<double>&,
+ double &RminLag, double &RminBez, int depth, int which) const;
+ void _computeRmin3(const fullMatrix<double>&, const fullVector<double>&,
+ double &RminLag, double &RminBez, int depth, bool debug = false) const;
+
+ double _subdivideForRmin(MetricData*, double RminLag, double tol, int which) const;
+
+ double _minp(const fullMatrix<double>&) const;
+ double _minp2(const fullMatrix<double>&) const;
+ double _minq(const fullMatrix<double>&) const;
+ double _maxp(const fullMatrix<double>&) const;
+ double _maxq(const fullMatrix<double>&) const;
+ void _minMaxA2(const fullMatrix<double>&, double &min, double &max) const;
+ void _minJ2P3(const fullMatrix<double>&, const fullVector<double>&, double &min) const;
+ void _maxAstK(const fullMatrix<double>&, const fullVector<double>&,
+ double minK, double a1, double &maxa) const; //wrong
+ void _maxAstK2(const fullMatrix<double>&, const fullVector<double>&,
+ double minK, double beta, double &maxa) const; //wrong
+ void _maxAstK3(const fullMatrix<double>&, const fullVector<double>&,
+ double minK, double beta, double &maxa) const; //poor bound
+ void _maxAstK4(const fullMatrix<double>&, const fullVector<double>&,
+ double minK, double beta, double &maxa) const; //better bound, WI positive ?
+ void _maxKstA(const fullMatrix<double>&, const fullVector<double>&,
+ double mina, double &maxK) const; //wrong
+ void _maxKstA2(const fullMatrix<double>&, const fullVector<double>&,
+ double mina, double beta, double &maxK) const; //faster
+ void _maxKstA3(const fullMatrix<double>&, const fullVector<double>&,
+ double mina, double beta, double &maxK) const; //better bound
+ void _minMaxJacobianSqr(const fullVector<double>&, double &min, double &max) const;
+
+private:
+ class gterIneq {
+ public:
+ bool operator()(const IneqData &id1, const IneqData &id2) const {
+ return id1.val > id2.val;
+ }
+ };
struct lessMinB {
bool operator()(const MetricData *md1, const MetricData *md2) const {
return md1->_RminBez > md2->_RminBez;
}
};
+};
- class IneqData {
+class MetricCoefficient {
+private:
+ class MetricData {
public:
- int i, j, k;
- double val;
+ fullMatrix<double> *_metcoeffs;
+ fullVector<double> *_jaccoeffs;
+ double _RminBez;
+ int _depth;
public:
- IneqData(double val, int i, int j, int k = -1) : i(i), j(j), k(k), val(val) {}
+ MetricData(fullMatrix<double> *m, fullVector<double> *j, double r, int d) :
+ _metcoeffs(m), _jaccoeffs(j), _RminBez(r), _depth(d) {}
+ ~MetricData() {
+ delete _metcoeffs;
+ delete _jaccoeffs;
+ }
};
- class gterIneq {
- public:
- bool operator()(const IneqData &id1, const IneqData &id2) const {
- return id1.val > id2.val;
+ struct lessMinB {
+ bool operator()(const MetricData *md1, const MetricData *md2) const {
+ return md1->_RminBez > md2->_RminBez;
}
};
@@ -54,9 +150,9 @@ private:
const GradientBasis *_gradients;
const bezierBasis *_bezier;
fullMatrix<double> _coefficientsLag, _coefficientsBez;
- std::map<int, std::vector<IneqData> > _ineqJ2, _ineqP3, _ineqA;
int __maxdepth, __numSubdivision;
std::vector<int> __numSub;
+ MetricBasis *_basis;
public:
MetricCoefficient(MElement*);
@@ -65,25 +161,9 @@ private:
void interpolate(const double *uvw, double *minmaxQ);
double getBoundRmin(double tol, int which);
- static int metricOrder(int tag);
-
private:
void _computeBezCoeff();
- void _fillInequalities(int order);
- void _lightenInequalities(double, int&, int&, int&); //TODO change
void _interpolateBezierPyramid(const double *uvw, double *minmaxQ);
- void _computeRmin1(fullMatrix<double>&, fullVector<double>&,
- double &RminLag, double &RminBez, int) const;
- void _computeRmin2(fullMatrix<double>&, fullVector<double>&,
- double &RminLag, double &RminBez, int depth, int which) const;
- double _subdivideForRmin(MetricData*, double RminLag, double tol, int which) const;
- double _minp(const fullMatrix<double>&) const;
- double _minq(const fullMatrix<double>&) const;
- double _maxp(const fullMatrix<double>&) const;
- double _maxq(const fullMatrix<double>&) const;
- void _minMaxA2(const fullMatrix<double>&, double &min, double &max) const;
- void _minMaxJacobianSqr(const fullVector<double>&, double &min, double &max) const;
- void _minJ2P3(const fullMatrix<double>&, const fullVector<double>&, double &min) const;
};
#endif