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

old developments

parent 75cd51b9
No related branches found
No related tags found
No related merge requests found
...@@ -202,12 +202,12 @@ double MElement::rhoShapeMeasure() ...@@ -202,12 +202,12 @@ double MElement::rhoShapeMeasure()
double MElement::metricShapeMeasure() double MElement::metricShapeMeasure()
{ {
return MetricBasis::minRCorner(this); return MetricBasis::getMinRCorner(this);
} }
double MElement::metricShapeMeasure2() double MElement::metricShapeMeasure2()
{ {
return MetricBasis::boundMinR(this); return MetricBasis::getMinR(this);
} }
double MElement::maxDistToStraight() const double MElement::maxDistToStraight() const
......
...@@ -562,3 +562,7 @@ int ElementType::getTag(int parentTag, int order, bool serendip) ...@@ -562,3 +562,7 @@ int ElementType::getTag(int parentTag, int order, bool serendip)
default : Msg::Warning("unknown element type %i, returning 0", parentTag); return 0; default : Msg::Warning("unknown element type %i, returning 0", parentTag); return 0;
} }
} }
int ElementType::getPrimaryTag(int tag) {
return getTag(ParentTypeFromTag(tag), 1);
}
...@@ -21,6 +21,9 @@ namespace ElementType ...@@ -21,6 +21,9 @@ namespace ElementType
// Give element tag from type, order & serendip // Give element tag from type, order & serendip
int getTag(int parentTag, int order, bool serendip = false); int getTag(int parentTag, int order, bool serendip = false);
// Give first order element tag
int getPrimaryTag(int tag);
} }
#endif #endif
...@@ -6,22 +6,40 @@ ...@@ -6,22 +6,40 @@
#include "MElement.h" #include "MElement.h"
#include "AnalyseCurvedMesh.h" #include "AnalyseCurvedMesh.h"
#include "MetricBasis.h" #include "MetricBasis.h"
#include "bezierBasis.h"
#include "JacobianBasis.h"
#include "BasisFactory.h" #include "BasisFactory.h"
#include "pointsGenerators.h" #include "pointsGenerators.h"
#include "MHexahedron.h"
#include "OS.h" #include "OS.h"
#include <queue> #include <queue>
#include <sstream> #include <sstream>
#include <limits> #include <limits>
#include <iomanip>
#include <cmath>
double MetricBasis::_tol = 1e-3; double MetricBasis::_tol = 1e-3;
//TODO Renommer fonctions plus explicitement (minmaxA,...) et rendre statique certaines fonctions
namespace { namespace {
double cubicCardanoRoot(double p, double q) double cubicCardanoRoot(double p, double q)
{ {
// solve the equation t^3 + p*t + q = 0
double A = q*q/4 + p*p*p/27; double A = q*q/4 + p*p*p/27;
if (A > 0) { if (A > 0) {
double sq = std::sqrt(A); double sq = std::sqrt(A);
return std::pow(-q/2+sq, 1/3.) + std::pow(-q/2-sq, 1/3.); double t1 = -q/2+sq, t2 = -q/2-sq;;
if (t1 < 0)
t1 = -std::pow(-t1, 1/3.);
else
t1 = std::pow(t1, 1/3.);
if (t2 < 0)
t2 = -std::pow(-t2, 1/3.);
else
t2 = std::pow(t2, 1/3.);
return t1 + t2;
} }
else { else {
double module = std::sqrt(-p*p*p/27); double module = std::sqrt(-p*p*p/27);
...@@ -58,6 +76,7 @@ MetricBasis::MetricBasis(int tag) : ...@@ -58,6 +76,7 @@ MetricBasis::MetricBasis(int tag) :
_jacobian(NULL), _type(ElementType::ParentTypeFromTag(tag)), _jacobian(NULL), _type(ElementType::ParentTypeFromTag(tag)),
_dim(ElementType::DimensionFromTag(tag)) _dim(ElementType::DimensionFromTag(tag))
{ {
//Msg::Fatal("Verifie si cette fonction est ok 1");
const bool serendip = false; const bool serendip = false;
const int metOrder = metricOrder(tag); const int metOrder = metricOrder(tag);
const int jacOrder = 3*metOrder/2; const int jacOrder = 3*metOrder/2;
...@@ -69,7 +88,6 @@ MetricBasis::MetricBasis(int tag) : ...@@ -69,7 +88,6 @@ MetricBasis::MetricBasis(int tag) :
else else
metricSpace = new FuncSpaceData(true, tag, false, metOrder+2, metricSpace = new FuncSpaceData(true, tag, false, metOrder+2,
metOrder, &serendip); metOrder, &serendip);
_gradients = BasisFactory::getGradientBasis(*metricSpace); _gradients = BasisFactory::getGradientBasis(*metricSpace);
_bezier = BasisFactory::getBezierBasis(*metricSpace); _bezier = BasisFactory::getBezierBasis(*metricSpace);
delete metricSpace; delete metricSpace;
...@@ -80,10 +98,15 @@ MetricBasis::MetricBasis(int tag) : ...@@ -80,10 +98,15 @@ MetricBasis::MetricBasis(int tag) :
jacSpace = new FuncSpaceData(true, tag, jacOrder, &serendip); jacSpace = new FuncSpaceData(true, tag, jacOrder, &serendip);
} }
else if (_type == TYPE_PYR) { else if (_type == TYPE_PYR) {
// The square of the jacobian must be in the same space than the cube of
// the metric, so +3
jacSpace = new FuncSpaceData(true, tag, false, jacOrder+3, jacSpace = new FuncSpaceData(true, tag, false, jacOrder+3,
jacOrder, &serendip); jacOrder, &serendip);
} }
else if (_type != TYPE_TRI && _type != TYPE_QUA) else if (_type == TYPE_TRI || _type == TYPE_QUA) {
// jacobian not needed
}
else
Msg::Fatal("metric not implemented for element tag %d", tag); Msg::Fatal("metric not implemented for element tag %d", tag);
if (jacSpace) { if (jacSpace) {
...@@ -94,59 +117,125 @@ MetricBasis::MetricBasis(int tag) : ...@@ -94,59 +117,125 @@ MetricBasis::MetricBasis(int tag) :
_fillInequalities(metOrder); _fillInequalities(metOrder);
} }
double MetricBasis::boundMinR(MElement *el) double MetricBasis::getMinR(MElement *el)
{ {
/*static int a = 0;
if (++a == 1) {
double pinf = 1/.0;
double minf = -1/.0;
double mynan = .0/.0;
Msg::Info("pinf=%g minf=%g mynan=%g", pinf, minf, mynan);
Msg::Info(". == . => %d %d %d", pinf == pinf, minf == minf, mynan == mynan);
Msg::Info("./. => %g %g %g", pinf/pinf, minf/minf, mynan/mynan);
Msg::Info("1. < . => %d %d %d", 1. < pinf, 1. < minf, 1. < mynan);
Msg::Info("1. > . => %d %d %d", 1. > pinf, 1. > minf, 1. > mynan);
Msg::Info("max => %g %g %g", std::max(1., pinf), std::max(1., minf), std::max(1., mynan));
Msg::Info("min => %g %g %g", std::min(1., pinf), std::min(1., minf), std::min(1., mynan));
Msg::Info("inv => %g %g %g", 1./pinf, 1./minf, 1./mynan);
Msg::Info("pinf>. => %d %d %d", pinf > pinf, pinf > minf, pinf > mynan);
Msg::Info(" ");
double giant = std::numeric_limits<double>::max();
double infinity = std::numeric_limits<double>::infinity();
Msg::Info("giant=%g, giant*10=%g infinity=%g", giant, giant*10, infinity);
Msg::Info("a=pinf, (a - 1) / (a + 1) => %g/%g => %g", pinf - 1, pinf + 1, (pinf - 1) / (pinf + 1));
Msg::Info("a=giant, (a - 1) / (a + 1) => %g/%g => %g", giant - 1, giant + 1, (giant - 1) / (giant + 1));
Msg::Info("giant*10 == infinity => %d", giant*10 == std::numeric_limits<double>::infinity());
Msg::Info("sqrt(pinf) => %g", std::sqrt(pinf));
Msg::Info("!(mynan >= 1) => %d, (mynan < 1) => %d", !(mynan >= 1), mynan < 1);
}*/
MetricBasis *metric = (MetricBasis*)BasisFactory::getMetricBasis(el->getTypeForMSH()); MetricBasis *metric = (MetricBasis*)BasisFactory::getMetricBasis(el->getTypeForMSH());
return metric->getBoundMinR(el); double m = metric->computeMinR(el);
//Msg::Info("returning %g", m);
return m;
} }
double MetricBasis::minRCorner(MElement *el) double MetricBasis::computeMinR(MElement *el) const
{ {
int tag = el->getTypeForMSH(); int nSampPnts = _gradients->getNumSamplingPoints();
int order = 1; int nMapping = _gradients->getNumMapNodes();
if (el->getType() == TYPE_TRI || el->getType() == TYPE_TET) order = 0;
const GradientBasis *gradients;
const JacobianBasis *jacobian;
if (el->getType() != TYPE_PYR) {
gradients = BasisFactory::getGradientBasis(tag, order);
jacobian = BasisFactory::getJacobianBasis(tag, order);
}
else {
const bool serendip = false;
FuncSpaceData fsd(true, tag, false, 1, 0, &serendip);
gradients = BasisFactory::getGradientBasis(fsd);
jacobian = BasisFactory::getJacobianBasis(fsd);
}
int nSampPnts = jacobian->getNumJacNodes();
if (el->getType() == TYPE_PYR) nSampPnts = 4;
int nMapping = gradients->getNumMapNodes();
// Nodes // Nodes
fullMatrix<double> nodes(nMapping, 3); fullMatrix<double> nodes(nMapping, 3);
el->getNodesCoord(nodes); el->getNodesCoord(nodes);
// Jacobian coefficients // Jacobian coefficients
fullVector<double> jacLag(jacobian->getNumJacNodes()); fullVector<double> *jac = NULL;
jacobian->getSignedJacobian(nodes, jacLag); if (_jacobian) {
fullVector<double> jacLag(_jacobian->getNumJacNodes());
jac = new fullVector<double>(_jacobian->getNumJacNodes());
_jacobian->getSignedIdealJacobian(nodes, jacLag);
_jacobian->lag2Bez(jacLag, *jac);
}
// Metric coefficients // Metric coefficients
fullMatrix<double> metCoeffLag; fullMatrix<double> metCoeffLag;
_fillCoeff<false>(el->getDim(), gradients, nodes, metCoeffLag); _fillCoeff<true>(el->getDim(), _gradients, nodes, metCoeffLag);
fullMatrix<double> *metCoeff;
metCoeff = new fullMatrix<double>(nSampPnts, metCoeffLag.size2());
_bezier->matrixLag2Bez.mult(metCoeffLag, *metCoeff);
// Compute min_corner(R) /*static double inf = std::numeric_limits<double>::infinity();
return _computeMinlagR(jacLag, metCoeffLag, nSampPnts); double minq = inf, maxq = -inf;
for (int k = 0; k < nSampPnts; ++k) {
minq = std::min(minq, (*metCoeff)(k, 0));
maxq = std::max(maxq, (*metCoeff)(k, 0));
}
double minp = 0, maxp = 0;
for (int i = 1; i < 7; ++i) {
double min = inf, max = -inf;
for (int k = 0; k < nSampPnts; ++k) {
min = std::min(min, (*metCoeff)(k, i));
max = std::max(max, (*metCoeff)(k, i));
}
if (min < 0 && max < 0) {
double tmp = min;
min = -max;
max = -tmp;
}
minp += min*min;
maxp += max*max;
}
minp = std::sqrt(minp);
maxp = std::sqrt(maxp);
double mina = minq/maxp;
double maxa = maxq/minp;
if (maxa-mina < _tol) {
return getMinSampledR(el, 0);
}*/
// Compute min(R, _tol)
double RminLag, RminBez;
//Msg::Info("element %d", el->getNum());
_computeBoundsRmin(*metCoeff, *jac, RminLag, RminBez);
//statsForMatlab(el, 20, NULL);
//Msg::Info("%.10g %.10g", RminBez, RminLag);
if (RminLag-RminBez < _tol) {
delete jac;
delete metCoeff;
return RminBez;
}
else {
//Msg::Info("Subdivision (%g vs %g)", RminLag, RminBez);
MetricData *md2 = new MetricData(metCoeff, jac, RminBez, 0, 0);
double minBez = _subdivideForRmin(md2, RminLag, _tol, el);
//Msg::Info("bez lag %.10g %.10g", minBez, RminLag);
return minBez;
}
} }
double MetricBasis::minSampledR(MElement *el, int order) double MetricBasis::getMinSampledR(MElement *el, int order)
{ {
MetricBasis *metric = (MetricBasis*)BasisFactory::getMetricBasis(el->getTypeForMSH()); MetricBasis *metric = (MetricBasis*)BasisFactory::getMetricBasis(el->getTypeForMSH());
return metric->getMinSampledR(el, order); return metric->computeMinSampledR(el, order);
} }
double MetricBasis::getMinSampledR(MElement *el, int deg) const double MetricBasis::computeMinSampledR(MElement *el, int deg) const
{ {
//Msg::Fatal("Verifie si cette fonction est ok b");
fullMatrix<double> samplingPoints; fullMatrix<double> samplingPoints;
bool serendip = false; bool serendip = false;
gmshGeneratePoints(FuncSpaceData(el, deg, &serendip), samplingPoints); gmshGeneratePoints(FuncSpaceData(el, deg, &serendip), samplingPoints);
...@@ -170,49 +259,118 @@ double MetricBasis::getMinSampledR(MElement *el, int deg) const ...@@ -170,49 +259,118 @@ double MetricBasis::getMinSampledR(MElement *el, int deg) const
return min; return min;
} }
double MetricBasis::getBoundMinR(MElement *el) const double MetricBasis::getMaxSampledR(MElement *el, int order)
{ {
int nSampPnts = _gradients->getNumSamplingPoints(); MetricBasis *metric = (MetricBasis*)BasisFactory::getMetricBasis(el->getTypeForMSH());
int nMapping = _gradients->getNumMapNodes(); return metric->computeMaxSampledR(el, order);
}
double MetricBasis::computeMaxSampledR(MElement *el, int deg) const
{
fullMatrix<double> samplingPoints;
bool serendip = false;
gmshGeneratePoints(FuncSpaceData(el, deg, &serendip), samplingPoints);
MetricData *md;
_getMetricData(el, md);
fullMatrix<double> R;
interpolate(el, md, samplingPoints, R);
if (R.size1() < 1) {
delete md;
return -1;
}
double max = R(0, 1);
for (int i = 1; i < R.size1(); ++i)
max = std::max(max, R(i, 1));
delete md;
return max;
}
void MetricBasis::lag2Bez(const fullMatrix<double> &metCoeffLag,
fullMatrix<double> &metCoeffBez) const
{
_bezier->matrixLag2Bez.mult(metCoeffLag, metCoeffBez);
}
double MetricBasis::getMinSampledGlobalRatio(MElement *el, int order)
{
MetricBasis *metric = (MetricBasis*)BasisFactory::getMetricBasis(el->getTypeForMSH());
return metric->computeMinSampledGlobalRatio(el, order);
}
double MetricBasis::computeMinSampledGlobalRatio(MElement *el, int deg) const
{
Msg::Fatal("Verifie si cette fonction est ok d");
fullMatrix<double> samplingPoints;
bool serendip = false;
gmshGeneratePoints(FuncSpaceData(el, deg, &serendip), samplingPoints);
MetricData *md;
_getMetricData(el, md);
fullMatrix<double> R;
interpolate(el, md, samplingPoints, R);
double min = R(0, 2);
double max = R(0, 3);
for (int i = 1; i < R.size1(); ++i) {
min = std::min(min, R(i, 2));
max = std::max(max, R(i, 3));
}
//Msg::Info("%g %g", min, max);
delete md;
return std::sqrt(min/max);
}
double MetricBasis::getMinRCorner(MElement *el)
{
Msg::Fatal("Verifie si cette fonction est ok e");
int tag = el->getTypeForMSH();
int order = 1;
if (el->getType() == TYPE_TRI || el->getType() == TYPE_TET) order = 0;
const GradientBasis *gradients;
const JacobianBasis *jacobian;
if (el->getType() != TYPE_PYR) {
gradients = BasisFactory::getGradientBasis(tag, order);
jacobian = BasisFactory::getJacobianBasis(tag, order);
}
else {
const bool serendip = false;
FuncSpaceData fsd(true, tag, false, 1, 0, &serendip);
gradients = BasisFactory::getGradientBasis(fsd);
jacobian = BasisFactory::getJacobianBasis(fsd);
}
int nSampPnts = jacobian->getNumJacNodes();
if (el->getType() == TYPE_PYR) nSampPnts = 4;
int nMapping = gradients->getNumMapNodes();
// Nodes // Nodes
fullMatrix<double> nodes(nMapping, 3); fullMatrix<double> nodes(nMapping, 3);
el->getNodesCoord(nodes); el->getNodesCoord(nodes);
// Jacobian coefficients // Jacobian coefficients
fullVector<double> *jac = NULL; fullVector<double> jacLag(jacobian->getNumJacNodes());
if (_jacobian) { jacobian->getSignedJacobian(nodes, jacLag);
fullVector<double> jacLag(_jacobian->getNumJacNodes());
jac = new fullVector<double>(_jacobian->getNumJacNodes());
_jacobian->getSignedIdealJacobian(nodes, jacLag);
_jacobian->lag2Bez(jacLag, *jac);
}
// Metric coefficients // Metric coefficients
fullMatrix<double> metCoeffLag; fullMatrix<double> metCoeffLag;
_fillCoeff<true>(el->getDim(), _gradients, nodes, metCoeffLag); _fillCoeff<false>(el->getDim(), gradients, nodes, metCoeffLag);
fullMatrix<double> *metCoeff;
metCoeff = new fullMatrix<double>(nSampPnts, metCoeffLag.size2());
_bezier->matrixLag2Bez.mult(metCoeffLag, *metCoeff);
// Compute min(R, _tol)
double RminLag, RminBez;
_computeRmin(*metCoeff, *jac, RminLag, RminBez);
if (RminLag-RminBez < MetricBasis::_tol) { // Compute min_corner(R)
delete jac; return _computeMinlagR(jacLag, metCoeffLag, nSampPnts);
delete metCoeff;
return RminBez;
}
else {
MetricData *md2 = new MetricData(metCoeff, jac, RminBez, 0, 0);
double tt = _subdivideForRmin(md2, RminLag, MetricBasis::_tol);
return tt;
}
} }
void MetricBasis::_fillInequalities(int metricOrder) void MetricBasis::_fillInequalities(int metricOrder)
{ {
//Msg::Fatal("Verifie si cette fonction est ok");
if (_type == TYPE_PYR) { if (_type == TYPE_PYR) {
_fillInequalitiesPyr(metricOrder); _fillInequalitiesPyr(metricOrder);
return; return;
...@@ -301,9 +459,6 @@ void MetricBasis::_fillInequalities(int metricOrder) ...@@ -301,9 +459,6 @@ void MetricBasis::_fillInequalities(int metricOrder)
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*metricOrder+1, l); hash += (exp(i, l)+exp(j, l)+exp(k, l)) * pow_int(3*metricOrder+1, l);
} }
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)); _ineqP3[hash].push_back(IneqData(num/den, i, j, k));
} }
} }
...@@ -361,6 +516,7 @@ void MetricBasis::_fillInequalities(int metricOrder) ...@@ -361,6 +516,7 @@ void MetricBasis::_fillInequalities(int metricOrder)
void MetricBasis::_fillInequalitiesPyr(int metricOrder) void MetricBasis::_fillInequalitiesPyr(int metricOrder)
{ {
//Msg::Fatal("Verifie si cette fonction est ok");
fullMatrix<int> exp(_bezier->_exponents.size1(), _bezier->_exponents.size2()); fullMatrix<int> exp(_bezier->_exponents.size1(), _bezier->_exponents.size2());
for (int i = 0; i < _bezier->_exponents.size1(); ++i) { for (int i = 0; i < _bezier->_exponents.size1(); ++i) {
for (int j = 0; j < _bezier->_exponents.size2(); ++j) { for (int j = 0; j < _bezier->_exponents.size2(); ++j) {
...@@ -424,9 +580,6 @@ void MetricBasis::_fillInequalitiesPyr(int metricOrder) ...@@ -424,9 +580,6 @@ void MetricBasis::_fillInequalitiesPyr(int metricOrder)
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*metricOrder+1, l); hash += (exp(i, l)+exp(j, l)+exp(k, l)) * pow_int(3*metricOrder+1, l);
} }
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)); _ineqP3[hash].push_back(IneqData(num/den, i, j, k));
} }
} }
...@@ -471,6 +624,7 @@ void MetricBasis::_fillInequalitiesPyr(int metricOrder) ...@@ -471,6 +624,7 @@ void MetricBasis::_fillInequalitiesPyr(int metricOrder)
void MetricBasis::_lightenInequalities(int &countj, int &countp, int &counta) void MetricBasis::_lightenInequalities(int &countj, int &countp, int &counta)
{ {
//Msg::Fatal("Verifie si cette fonction est ok");
double tol = .0; double tol = .0;
std::map<int, std::vector<IneqData> >::iterator it, itbeg[3], itend[3]; std::map<int, std::vector<IneqData> >::iterator it, itbeg[3], itend[3];
...@@ -506,10 +660,74 @@ void MetricBasis::_lightenInequalities(int &countj, int &countp, int &counta) ...@@ -506,10 +660,74 @@ void MetricBasis::_lightenInequalities(int &countj, int &countp, int &counta)
bool MetricBasis::validateBezierForMetricAndJacobian() bool MetricBasis::validateBezierForMetricAndJacobian()
{ {
Msg::Fatal("Verifie si cette fonction est ok f");
Msg::Info("Testing Bezier interpolation and subdivision " Msg::Info("Testing Bezier interpolation and subdivision "
"for jacobien and metric on all element types..."); "for jacobien and metric on all element types...");
int numError = 0; int numError = 0;
// Parameters
const int numElem = 10000000; ////////
MElement **elements;
elements = new MElement*[numElem];
for (int tag = 1; tag <= MSH_NUM_TYPE; ++tag) {
if (tag != 9) continue; ////////
const int numPrimNode = 3; ////////
const int order = ElementType::OrderFromTag(tag);
const int dim = ElementType::DimensionFromTag(tag);
const int primTag = ElementType::getPrimaryTag(tag);
Msg::Info("... testing element tag %d", tag);
// Get reference nodes
const nodalBasis *mapBasis = BasisFactory::getNodalBasis(tag);
fullMatrix<double> nodes;
mapBasis->getReferenceNodes(nodes);
// Create 'numElem' elements more and more randomized
for (int iel = 0; iel < numElem; ++iel) {
const double range = static_cast<double>(iel) / (numElem-1) / order;
const double rangePrim = static_cast<double>(iel/2000) / ((numElem/2000)-1) / order; ////////
const double rangeSub = static_cast<double>(iel%2000) / 1999 / order; ////////
if (!(iel%200)) Msg::Info("%g %g", rangePrim, rangeSub);
std::vector<MVertex*> vertices(nodes.size1());
for (int i = 0; i < numPrimNode; ++i) {
vertices[i] = new MVertex(nodes(i, 0) + symRand(rangePrim),
nodes(i, 1) + symRand(rangePrim),
nodes(i, 2) + symRand(rangePrim));
}
MElement *primEl = MElement::createElement(primTag, vertices);
for (int i = numPrimNode; i < nodes.size1(); ++i) {
SPoint3 p;
primEl->pnt(nodes(i, 0), nodes(i, 1), nodes(i, 2), p);
vertices[i] = new MVertex(p.x() + symRand(rangeSub),
p.y() + symRand(rangeSub),
p.z() + symRand(rangeSub));
}
MElement *el = MElement::createElement(tag, vertices);
if (!el) {
Msg::Error("MElement was unable to create element for tag %d", tag);
++numError;
}
elements[iel] = el;
}
GMSH_AnalyseCurvedMeshPlugin plugin = GMSH_AnalyseCurvedMeshPlugin();
plugin.test(elements, numElem, dim);
}
//=================
return 987;
//=================
/*Msg::Info("Testing Bezier interpolation and subdivision "
"for jacobien and metric on all element types...");
int numError = 0;
// Parameters // Parameters
const int numType = 6; const int numType = 6;
const int acceptedTypes[numType] = {TYPE_TRI, TYPE_QUA, TYPE_TET, const int acceptedTypes[numType] = {TYPE_TRI, TYPE_QUA, TYPE_TET,
...@@ -566,7 +784,7 @@ bool MetricBasis::validateBezierForMetricAndJacobian() ...@@ -566,7 +784,7 @@ bool MetricBasis::validateBezierForMetricAndJacobian()
dim > 2 ? nodes(i, 2) + symRand(range) : 0); dim > 2 ? nodes(i, 2) + symRand(range) : 0);
} }
MElement *el = MElement::createElement(tag, vertices); MElement *el = MElement::createElement(tag, vertices);
if (el != NULL) { if (!el) {
Msg::Error("MElement was unable to create element for tag %d", tag); Msg::Error("MElement was unable to create element for tag %d", tag);
++numError; ++numError;
} }
...@@ -578,7 +796,7 @@ bool MetricBasis::validateBezierForMetricAndJacobian() ...@@ -578,7 +796,7 @@ bool MetricBasis::validateBezierForMetricAndJacobian()
if (numError) Msg::Error("Validation of Bezier terminated with %d errors!", numError); if (numError) Msg::Error("Validation of Bezier terminated with %d errors!", numError);
else Msg::Info("Validation of Bezier terminated without errors"); else Msg::Info("Validation of Bezier terminated without errors");
return numError; return numError;*/
} }
int MetricBasis::validateBezierForMetricAndJacobian(MElement *el, int MetricBasis::validateBezierForMetricAndJacobian(MElement *el,
...@@ -587,6 +805,7 @@ int MetricBasis::validateBezierForMetricAndJacobian(MElement *el, ...@@ -587,6 +805,7 @@ int MetricBasis::validateBezierForMetricAndJacobian(MElement *el,
double toleranceTensor, double toleranceTensor,
double tolerance) double tolerance)
{ {
Msg::Fatal("Verifie si cette fonction est ok g");
const int tag = el->getTypeForMSH(); const int tag = el->getTypeForMSH();
const MetricBasis *metricBasis = BasisFactory::getMetricBasis(tag); const MetricBasis *metricBasis = BasisFactory::getMetricBasis(tag);
...@@ -638,12 +857,13 @@ void MetricBasis::interpolate(const MElement *el, ...@@ -638,12 +857,13 @@ void MetricBasis::interpolate(const MElement *el,
const fullMatrix<double> &nodes, const fullMatrix<double> &nodes,
fullMatrix<double> &R) const fullMatrix<double> &R) const
{ {
//Msg::Fatal("Verifie si cette fonction est ok h");
fullMatrix<double> &metcoeffs = *md->_metcoeffs, *metric = new fullMatrix<double>; fullMatrix<double> &metcoeffs = *md->_metcoeffs, *metric = new fullMatrix<double>;
fullVector<double> &jaccoeffs = *md->_jaccoeffs, *jac = new fullVector<double>; fullVector<double> &jaccoeffs = *md->_jaccoeffs, *jac = new fullVector<double>;
_bezier->interpolate(metcoeffs, nodes, *metric); _bezier->interpolate(metcoeffs, nodes, *metric);
R.resize(nodes.size1(), 2); R.resize(nodes.size1(), 4);
switch (metcoeffs.size2()) { switch (metcoeffs.size2()) {
case 1: case 1:
...@@ -652,6 +872,8 @@ void MetricBasis::interpolate(const MElement *el, ...@@ -652,6 +872,8 @@ void MetricBasis::interpolate(const MElement *el,
break; break;
case 3: case 3:
((MetricBasis*)this)->file << "q p myR eigR" << std::endl;
((MetricBasis*)this)->file << _bezier->getNumLagCoeff() << std::endl;
for (int i = 0; i < R.size1(); ++i) { for (int i = 0; i < R.size1(); ++i) {
// Compute from q, p // Compute from q, p
double p = pow((*metric)(i, 1), 2); double p = pow((*metric)(i, 1), 2);
...@@ -667,10 +889,18 @@ void MetricBasis::interpolate(const MElement *el, ...@@ -667,10 +889,18 @@ void MetricBasis::interpolate(const MElement *el,
fullMatrix<double> vecLeft(2, 2), vecRight(2, 2); fullMatrix<double> vecLeft(2, 2), vecRight(2, 2);
metricTensor.eig(valReal, valImag, vecLeft, vecRight, true); metricTensor.eig(valReal, valImag, vecLeft, vecRight, true);
R(i, 1) = std::sqrt(valReal(0) / valReal(1)); R(i, 1) = std::sqrt(valReal(0) / valReal(1));
R(i, 2) = valReal(0);
R(i, 3) = valReal(1);
((MetricBasis*)this)->file << (*metric)(i, 0) << " ";
((MetricBasis*)this)->file << p << " ";
((MetricBasis*)this)->file << R(i, 0) << " ";
((MetricBasis*)this)->file << R(i, 1) << std::endl;
} }
break; break;
case 7: case 7:
((MetricBasis*)this)->file << "a K myR eigR" << std::endl;
((MetricBasis*)this)->file << _bezier->getNumLagCoeff() << std::endl;
_jacobian->getBezier()->interpolate(jaccoeffs, nodes, *jac); _jacobian->getBezier()->interpolate(jaccoeffs, nodes, *jac);
for (int i = 0; i < R.size1(); ++i) { for (int i = 0; i < R.size1(); ++i) {
// Compute from q, p, J // Compute from q, p, J
...@@ -692,6 +922,12 @@ void MetricBasis::interpolate(const MElement *el, ...@@ -692,6 +922,12 @@ void MetricBasis::interpolate(const MElement *el,
fullMatrix<double> vecLeft(3, 3), vecRight(3, 3); fullMatrix<double> vecLeft(3, 3), vecRight(3, 3);
metricTensor.eig(valReal, valImag, vecLeft, vecRight, true); metricTensor.eig(valReal, valImag, vecLeft, vecRight, true);
R(i, 1) = std::sqrt(valReal(0) / valReal(2)); R(i, 1) = std::sqrt(valReal(0) / valReal(2));
R(i, 2) = valReal(0);
R(i, 3) = valReal(2);
((MetricBasis*)this)->file << (*metric)(i, 0)/p << " ";
((MetricBasis*)this)->file << (*jac)(i)*(*jac)(i)/p/p/p << " ";
((MetricBasis*)this)->file << R(i, 0) << " ";
((MetricBasis*)this)->file << R(i, 1) << std::endl;
} }
break; break;
...@@ -710,6 +946,7 @@ void MetricBasis::interpolateAfterNSubdivisions( ...@@ -710,6 +946,7 @@ void MetricBasis::interpolateAfterNSubdivisions(
fullMatrix<double> &uvw, fullMatrix<double> &uvw,
fullMatrix<double> &metric) const fullMatrix<double> &metric) const
{ {
Msg::Fatal("Verifie si cette fonction est ok i");
// Interpolate metric after 'numSub' random subdivision at // Interpolate metric after 'numSub' random subdivision at
// 'numPnt' random points. // 'numPnt' random points.
// Return: isub, the subdomain tag of each subdivision, // Return: isub, the subdomain tag of each subdivision,
...@@ -798,6 +1035,8 @@ void MetricBasis::interpolateAfterNSubdivisions( ...@@ -798,6 +1035,8 @@ void MetricBasis::interpolateAfterNSubdivisions(
void MetricBasis::statsForMatlab(MElement *el, int deg, MetricData *md) const void MetricBasis::statsForMatlab(MElement *el, int deg, MetricData *md) const
{ {
//return;
//Msg::Fatal("Verifie si cette fonction est ok h");
fullMatrix<double> samplingPoints; fullMatrix<double> samplingPoints;
bool serendip = false; bool serendip = false;
gmshGeneratePoints(FuncSpaceData(el, deg, &serendip), samplingPoints); gmshGeneratePoints(FuncSpaceData(el, deg, &serendip), samplingPoints);
...@@ -805,7 +1044,7 @@ void MetricBasis::statsForMatlab(MElement *el, int deg, MetricData *md) const ...@@ -805,7 +1044,7 @@ void MetricBasis::statsForMatlab(MElement *el, int deg, MetricData *md) const
if (!md) _getMetricData(el, md); if (!md) _getMetricData(el, md);
static unsigned int aa = 0; static unsigned int aa = 0;
if (md->_num < 100 && ++aa < 200) { if (md->_num < 100000 && ++aa < 1000) {
std::stringstream name; std::stringstream name;
name << "metricStat_" << el->getNum() << "_"; name << "metricStat_" << el->getNum() << "_";
name << (md->_num % 10); name << (md->_num % 10);
...@@ -817,42 +1056,47 @@ void MetricBasis::statsForMatlab(MElement *el, int deg, MetricData *md) const ...@@ -817,42 +1056,47 @@ void MetricBasis::statsForMatlab(MElement *el, int deg, MetricData *md) const
((MetricBasis*)this)->file.open(name.str().c_str(), std::fstream::out); ((MetricBasis*)this)->file.open(name.str().c_str(), std::fstream::out);
{ {
int dim = el->getDim(); ((MetricBasis*)this)->file << "mina minKfast minKsharp cmin beta_m beta_M beta c_m" << std::endl;
double mina, maxa, minKs, minKf, c, beta_m, beta_M, beta, c_m;
fullMatrix<double> *coeff = md->_metcoeffs; fullMatrix<double> *coeff = md->_metcoeffs;
fullVector<double> *jac = md->_jaccoeffs; fullVector<double> *jac = md->_jaccoeffs;
double mina, maxa, minK;
_minMaxA(*coeff, mina, maxa); _minMaxA(*coeff, mina, maxa);
double dRda, term1, phip, beta; if (el->getDim() == 3) {
double maxaPos, maxaNeg; _minKfast(*coeff, *jac, minKf);
double maxKfast, maxKsharp; _minKsharp(*coeff, *jac, minKs);
if (dim == 3) { double minK = minKf;
_minK(*coeff, *jac, minK); _moveInsideDomain(mina, minK, true);
double tmpa = mina, tmpK = minK; double p_dRda, p_dRdK;
_computeTermBeta(tmpa, tmpK, dRda, term1, phip); _computePseudoDerivatives(mina, minK, p_dRda, p_dRdK);
beta = -3 * mina*mina * term1 / dRda / 6; double slope = -p_dRda/p_dRdK;
c = .5*(3*minK/mina - slope);
_maxAstKneg(*coeff, *jac, minK, beta, maxaPos); _computeBoundingCurve(*coeff, *jac, beta_m, c, true);
_maxAstKpos(*coeff, *jac, minK, beta, maxaNeg); _computeBoundingCurve(*coeff, *jac, beta_M, c, false);
beta = -p_dRda/(p_dRdK*3*mina*mina);
_maxKstAfast(*coeff, *jac, mina, beta, maxKfast); if (beta < 0)
_maxKstAsharp(*coeff, *jac, mina, beta, maxKsharp); _computeLowerBoundC(*coeff, *jac, beta, c_m);
else
c_m = -1;
} }
else { else {
minK = dRda = term1 = phip = beta = maxaPos = maxaNeg = maxKfast = maxKsharp = -1; minKs = minKf = c = beta_m = beta_M = beta = c_m = -1;
} }
((MetricBasis*)this)->file << mina << " " << maxa << " " << minK << " "; ((MetricBasis*)this)->file << std::setprecision(15);
((MetricBasis*)this)->file << dRda << " " << term1 << " " << phip << " "; ((MetricBasis*)this)->file << mina << " " << minKf << " " << minKs << " ";
((MetricBasis*)this)->file << beta << " "; ((MetricBasis*)this)->file << c << " " << beta_m << " " << beta_M << " ";
((MetricBasis*)this)->file << maxaPos << " " << maxaNeg << " "; ((MetricBasis*)this)->file << beta << " " << c_m << std::endl;
((MetricBasis*)this)->file << maxKfast << " " << maxKsharp << std::endl; ((MetricBasis*)this)->file << std::setprecision(8);
} }
} }
fullMatrix<double> R; fullMatrix<double> R;
interpolate(el, md, samplingPoints, R); interpolate(el, md, samplingPoints, R);
((MetricBasis*)this)->file.close();
return;
if (R.size1() < 1) { if (R.size1() < 1) {
((MetricBasis*)this)->file << -1 << " "; ((MetricBasis*)this)->file << -1 << " ";
...@@ -893,139 +1137,124 @@ int MetricBasis::metricOrder(int tag) ...@@ -893,139 +1137,124 @@ int MetricBasis::metricOrder(int tag)
} }
} }
void MetricBasis::_computeRmin( void MetricBasis::_computeBoundsRmin(
const fullMatrix<double> &coeff, const fullVector<double> &jac, const fullMatrix<double> &coeff, const fullVector<double> &jac,
double &RminLag, double &RminBez) const double &sqrtRminLag, double &sqrtRminBez) const
{ {
RminLag = _computeMinlagR(jac, coeff, _bezier->getNumLagCoeff()); sqrtRminLag = _computeMinlagR(jac, coeff, _bezier->getNumLagCoeff());
if (RminLag == 0) { if (sqrtRminLag <= _tol) {
RminBez = 0; sqrtRminBez = .0;
return; return;
} }
if (coeff.size2() == 3) { // 2d element if (coeff.size2() == 3) { // 2d element
double mina, dummy; double mina, dummy;
_minMaxA(coeff, mina, dummy); _minMaxA(coeff, mina, dummy);
RminBez = std::sqrt(_R2Dsafe(mina)); sqrtRminBez = std::sqrt(_R2Dsafe(mina));
return;
}
double minK;
_minK(coeff, jac, minK);
if (minK < 1e-10) {
RminBez = 0;
return; return;
} }
double mina, dummy; double minK, mina, dummy;
_minKfast(coeff, jac, minK);
_minMaxA(coeff, mina, dummy); _minMaxA(coeff, mina, dummy);
double term1, dRda, phip; _moveInsideDomain(mina, minK, true);
_computeTermBeta(mina, minK, dRda, term1, phip);
double p_dRda, p_dRdK;
_computePseudoDerivatives(mina, minK, p_dRda, p_dRdK);
if (dRda < 0) { if (p_dRda < 0) {
// TODO : better am ? //Msg::Info("case 1 & 2");
double amApprox, da; // cases 1 & 2 => compute a lower bounding curve K = beta * a^3 + c * a
// with c such that the curve that pass through (minK, mina) has the
// slope equal to -dRda/dRdK
double slope = -p_dRda/p_dRdK;
double c = .5*(3*minK/mina - slope);
double beta;
_computeBoundingCurve(coeff, jac, beta, c, true);
{ {
const double p = -3; double beta = (minK/mina-c)/mina/mina;
double q = -minK - 2; if (std::abs(p_dRda + p_dRdK * (3*beta*mina*mina+c)) > 1e-12) {
const double a1 = cubicCardanoRoot(p, q); Msg::Error("Derivative along curve not zero %g", p_dRda + p_dRdK * (3*beta*mina*mina+c));
const double phim = std::acos(-1/a1) - M_PI/3; }
q = -minK + 2*std::cos(3*phim);
amApprox = cubicCardanoRoot(p, q);
if (minK < 10)
da = -.3;
else if (minK < 20)
da = -.25;
else if (minK < 35)
da = -.2;
else if (minK < 70)
da = -.15;
else if (minK < 175)
da = -.1;
else
da = -.05;
} }
double beta = -3 * mina*mina * term1 / dRda / 6; double a_int = mina, K_int = minK;
double maxa; bool onVertical = _intersectionCurveLeftCorner(beta, c, a_int, K_int);
if (beta*minK-mina*mina*mina < 0) if (onVertical) {
_maxAstKneg(coeff, jac, minK, beta, maxa); // then the minimum is on the curve, we prefer to return this
else { sqrtRminBez = std::sqrt(_R3Dsafe(mina, minK));
_maxAstKpos(coeff, jac, minK, beta, maxa); //Msg::Info("onvertical");
if (maxa < amApprox && beta*minK-maxa*maxa*maxa < 0) return;
_maxAstKneg(coeff, jac, minK, beta, maxa);
} }
/*if (a_int - mina < ???) {
sqrtRminBez = std::sqrt(_R3Dsafe(a_int, minK));
return;
}*/
maxa = std::max(mina, maxa); bool haveToCompute_a0;
if (amApprox*amApprox*amApprox+da < maxa*maxa*maxa) { if (_moveInsideDomain(a_int, K_int, true)) {
// compute better am haveToCompute_a0 = true;
//
double am0 = std::pow(amApprox*amApprox*amApprox+da, 1/3.);
double am1 = std::pow(amApprox*amApprox*amApprox+da+.05, 1/3.);
//double am0S = am0, am1S = am1;
double am = (am0 + am1)/2;
double R0 = _R3Dsafe(am0, minK);
double R1 = _R3Dsafe(am1, minK);
double Rnew = _R3Dsafe(am, minK);
int cnt = 0;
while (std::abs(R0-Rnew) > _tol*.01 || std::abs(R1-Rnew) > _tol*.01) {
++cnt;
if (R0 > R1) {
am0 = am;
R0 = Rnew;
} }
else { else {
am1 = am; double p_dRda, p_dRdK;
R1 = Rnew; _computePseudoDerivatives(a_int, K_int, p_dRda, p_dRdK);
if (p_dRda + p_dRdK * (3*beta*a_int*a_int+c) < -1e-5) {
Msg::Error("Derivative along curve not positive or zero %g", p_dRda + p_dRdK * (3*beta*a_int*a_int+c));
Msg::Info("(%g, %g) vs (%g, %g), diff (%g, %g)", mina, minK, a_int, K_int, a_int-mina, K_int-minK);
double beta2 = (minK/mina-c)/mina/mina;
Msg::Info("beta %g - %g = %g", beta, beta2, beta-beta2);
} }
am = (am0 + am1)/2; haveToCompute_a0 = p_dRda > 0;
Rnew = _R3Dsafe(am, minK);
} }
if (am < maxa) { if (haveToCompute_a0) {
RminBez = _R3Dsafe(am, minK); double a0 = _computeAbscissaMinR(mina, minK);
RminBez = std::sqrt(RminBez); sqrtRminBez = std::sqrt(_R3Dsafe(a0, minK));
return;
} }
else {
sqrtRminBez = std::sqrt(_R3Dsafe(a_int, K_int));
} }
RminBez = _R3Dsafe(maxa, minK);
RminBez = std::sqrt(RminBez);
return; return;
} }
else if (term1 < 0) { else if (p_dRdK < 0) {
double maxK; //Msg::Info("case 4 & 5");
double beta = -3 * mina*mina * term1 / dRda / 6; double slope = -p_dRda/p_dRdK;
if (beta*minK-mina*mina*mina > 0) Msg::Fatal("Arf pas prevu"); double c = .5*(3*minK/mina - slope);
//_maxKstAsharp(coeff, jac, mina, beta, maxK); double beta;
_maxKstAfast(coeff, jac, mina, beta, maxK); _computeBoundingCurve(coeff, jac, beta, c, false);
const double x = .5*(maxK-mina*mina*mina+3*mina); double a_int = mina, K_int = minK;
const double phimin = std::acos(-1/mina) - M_PI/3;
double myphi; if (_intersectionCurveLeftCorner(beta, c, a_int, K_int)) {
if (std::abs(x) > 1) { /*if (std::abs(a_int - mina) < _tol / 100 &&
myphi = phimin; std::abs(K_int - minK) < _tol / 100) {
sqrtRminBez = std::sqrt(_R3Dsafe(mina, minK));
return;
}*/
// We compute K0 and return _R3Dsafe(mina, min(K0,K_int))
double K0 = 2*std::cos(3*std::acos(-1/mina)-M_PI) + (mina*mina-3) * mina;
sqrtRminBez = std::sqrt(_R3Dsafe(mina, std::min(K0, K_int)));
} }
else { else {
const double phimaxK = std::acos(x)/3; if (_moveInsideDomain(a_int, K_int, false))
myphi = std::max(phimin, phimaxK); sqrtRminBez = std::sqrt(_R3Dsafe(a_int, K_int));
else
sqrtRminBez = std::sqrt(_computeMinRAlongCurve(beta, c, a_int, -1));
} }
RminBez = (mina+2*std::cos(myphi+2*M_PI/3))/(mina+2*std::cos(myphi));
RminBez = std::sqrt(RminBez);
return; return;
} }
else { else {
RminBez = (mina+2*std::cos(phip+M_PI/3))/(mina+2*std::cos(phip-M_PI/3)); //Msg::Info("case 3");
RminBez = std::sqrt(RminBez); sqrtRminBez = std::sqrt(_R3Dsafe(mina, minK));
return; return;
} }
} }
void MetricBasis::_computeRmax( void MetricBasis::_computeBoundsRmax(
const fullMatrix<double> &coeff, const fullVector<double> &jac, const fullMatrix<double> &coeff, const fullVector<double> &jac,
double &RmaxLag) const double &RmaxLag) const
{ {
Msg::Fatal("Verifie si cette fonction est ok k");
RmaxLag = 0.; RmaxLag = 0.;
for (int i = 0; i < _bezier->getNumLagCoeff(); ++i) { for (int i = 0; i < _bezier->getNumLagCoeff(); ++i) {
...@@ -1046,8 +1275,9 @@ void MetricBasis::_computeRmax( ...@@ -1046,8 +1275,9 @@ void MetricBasis::_computeRmax(
} }
} }
double MetricBasis::_subdivideForRmin(MetricData *md, double RminLag, double tol, MElement *el) const double MetricBasis::_subdivideForRmin(MetricData *md, double &RminLag, double tol, MElement *el) const
{ {
//Msg::Fatal("Verifie si cette fonction est ok l");
std::priority_queue<MetricData*, std::vector<MetricData*>, lessMinB> subdomains; std::priority_queue<MetricData*, std::vector<MetricData*>, lessMinB> subdomains;
const bool for3d = md->_jaccoeffs; const bool for3d = md->_jaccoeffs;
const int numCoeff = md->_metcoeffs->size2(); const int numCoeff = md->_metcoeffs->size2();
...@@ -1058,6 +1288,7 @@ double MetricBasis::_subdivideForRmin(MetricData *md, double RminLag, double tol ...@@ -1058,6 +1288,7 @@ double MetricBasis::_subdivideForRmin(MetricData *md, double RminLag, double tol
std::vector<fullVector<double>*> trash; std::vector<fullVector<double>*> trash;
int k = 0;
while (RminLag - subdomains.top()->_RminBez > tol && subdomains.size() < 25000) { while (RminLag - subdomains.top()->_RminBez > tol && subdomains.size() < 25000) {
fullMatrix<double> *subcoeffs, *coeff; fullMatrix<double> *subcoeffs, *coeff;
fullVector<double> *subjac, *jac = NULL; fullVector<double> *subjac, *jac = NULL;
...@@ -1082,17 +1313,22 @@ double MetricBasis::_subdivideForRmin(MetricData *md, double RminLag, double tol ...@@ -1082,17 +1313,22 @@ double MetricBasis::_subdivideForRmin(MetricData *md, double RminLag, double tol
jac->setAsProxy(*subjac, i * numJacPnts, numJacPnts); jac->setAsProxy(*subjac, i * numJacPnts, numJacPnts);
} }
double minLag, minBez; double minLag, minBez;
_computeRmin(*coeff, *jac, minLag, minBez); _computeBoundsRmin(*coeff, *jac, minLag, minBez);
RminLag = std::min(RminLag, minLag); RminLag = std::min(RminLag, minLag);
int newNum = num + (i+1) * pow_int(10, depth); int newNum = num + (i+1) * pow_int(10, depth);
if (depth > 8) newNum = 999999999;
MetricData *metData = new MetricData(coeff, jac, minBez, depth+1, newNum); MetricData *metData = new MetricData(coeff, jac, minBez, depth+1, newNum);
//statsForMatlab(el, 15, metData);
//if (el) statsForMatlab(el, 20, metData); //if (el) statsForMatlab(el, 20, metData);
subdomains.push(metData); subdomains.push(metData);
} }
if (for3d) trash.push_back(subjac); if (for3d) trash.push_back(subjac);
delete subcoeffs; delete subcoeffs;
++k;
} }
//Msg::Info("num subdiv %d", k);
double ans = subdomains.top()->_RminBez; double ans = subdomains.top()->_RminBez;
while (subdomains.size()) { while (subdomains.size()) {
md = subdomains.top(); md = subdomains.top();
...@@ -1138,6 +1374,7 @@ template<bool ideal> ...@@ -1138,6 +1374,7 @@ template<bool ideal>
void MetricBasis::_fillCoeff(int dim, const GradientBasis *gradients, void MetricBasis::_fillCoeff(int dim, const GradientBasis *gradients,
const fullMatrix<double> &nodes, fullMatrix<double> &coeff) const fullMatrix<double> &nodes, fullMatrix<double> &coeff)
{ {
//Msg::Fatal("Verifie si cette fonction est ok o");
const int nSampPnts = gradients->getNumSamplingPoints(); const int nSampPnts = gradients->getNumSamplingPoints();
switch (dim) { switch (dim) {
...@@ -1159,15 +1396,11 @@ void MetricBasis::_fillCoeff(int dim, const GradientBasis *gradients, ...@@ -1159,15 +1396,11 @@ void MetricBasis::_fillCoeff(int dim, const GradientBasis *gradients,
for (int i = 0; i < nSampPnts; i++) { for (int i = 0; i < nSampPnts; i++) {
const double &dxdX = dxydX(i,0), &dydX = dxydX(i,1); const double &dxdX = dxydX(i,0), &dydX = dxydX(i,1);
const double &dxdY = dxydY(i,0), &dydY = dxydY(i,1); const double &dxdY = dxydY(i,0), &dydY = dxydY(i,1);
double dzdX, dzdY; double dzdX = 0, dzdY = 0;
if (nodes.size2() > 2) { if (nodes.size2() > 2) {
dzdX = dxydX(i,2); dzdX = dxydX(i,2);
dzdY = dxydY(i,2); dzdY = dxydY(i,2);
} }
else {
dzdX = 0;
dzdY = 0;
}
const double dvxdX = dxdX*dxdX + dydX*dydX + dzdX*dzdX; const double dvxdX = dxdX*dxdX + dydX*dydX + dzdX*dzdX;
const double dvxdY = dxdY*dxdY + dydY*dydY + dzdY*dzdY; const double dvxdY = dxdY*dxdY + dydY*dydY + dzdY*dzdY;
coeff(i, 0) = (dvxdX + dvxdY) / 2; coeff(i, 0) = (dvxdX + dvxdY) / 2;
...@@ -1218,46 +1451,43 @@ double MetricBasis::_computeMinlagR(const fullVector<double> &jac, ...@@ -1218,46 +1451,43 @@ double MetricBasis::_computeMinlagR(const fullVector<double> &jac,
for (int i = 0; i < num; ++i) { for (int i = 0; i < num; ++i) {
if (jac(i) <= 0.) return 0; if (jac(i) <= 0.) return 0;
const double q = coeff(i, 0); const double &q = coeff(i, 0);
double p = 0; double p = 0;
for (int k = 1; k < 7; ++k) { for (int k = 1; k < 7; ++k) {
p += pow_int(coeff(i, k), 2); p += pow_int(coeff(i, k), 2);
} }
p = std::sqrt(p); p = std::sqrt(p);
const double a = q/p;
if (a > 1e4) { // TODO: from _tol ? Rmin = std::min(Rmin, _R3Dsafe(q, p, jac(i)));
Rmin = std::min(Rmin, std::sqrt((a - std::sqrt(3.)) / (a + std::sqrt(3.))));
}
else {
const double tmpR = _R3Dsafe(a, jac(i)/p/p*jac(i)/p);
Rmin = std::min(Rmin, std::sqrt(tmpR));
}
} }
return Rmin; break;
case 3: case 3:
for (int i = 0; i < num; ++i) { for (int i = 0; i < num; ++i) {
const double &q = coeff(i, 0); const double &q = coeff(i, 0);
const double p = pow_int(coeff(i, 1), 2) + pow_int(coeff(i, 2), 2); double p = pow_int(coeff(i, 1), 2) + pow_int(coeff(i, 2), 2);
const double tmpR = _R2Dsafe(q, std::sqrt(p)); p = std::sqrt(p);
Rmin = std::min(Rmin, std::sqrt(tmpR));
Rmin = std::min(Rmin, _R2Dsafe(q, p));
} }
return Rmin; break;
default: default:
Msg::Error("coeff have not right number of column"); Msg::Error("coeff have not right number of column");
return -1; return -1;
} }
return std::sqrt(Rmin);
} }
void MetricBasis::_minMaxA( void MetricBasis::_minMaxA(const fullMatrix<double> &coeff,
const fullMatrix<double> &coeff, double &min, double &max) const double &mina, double &max) const
{ {
min = std::numeric_limits<double>::max(); // TODO calculer le max ?
min = 1e10; double upperBound = std::numeric_limits<double>::infinity();
//max = 1; // max is not used actually double lowerBound = -upperBound;
double minLowBound = -min;
std::map<int, std::vector<IneqData> >::const_iterator it = _ineqA.begin(); std::map<int, std::vector<IneqData> >::const_iterator it = _ineqA.begin();
while (it != _ineqA.end()) { while (it != _ineqA.end()) {
double num = 0; double num = 0;
double den = 0; double den = 0;
...@@ -1271,36 +1501,34 @@ void MetricBasis::_minMaxA( ...@@ -1271,36 +1501,34 @@ void MetricBasis::_minMaxA(
den += it->second[k].val * tmp; den += it->second[k].val * tmp;
num += it->second[k].val * coeff(i, 0) * coeff(j, 0); num += it->second[k].val * coeff(i, 0) * coeff(j, 0);
} }
double val = num/den;
if (num < 0) { // mina has to satisfy: mina^2 * den <= num
if (den > 0) { if (den == 0) {
_minA(coeff, min); if (num >= 0) {} // ok, nothing to do
else {// impossible to satisfy.
_minAfast(coeff, mina);
return; return;
} }
minLowBound = std::max(val, minLowBound);
}
else if (den > 0) {
min = std::min(val, min);
//max = std::max(val, max);
}
++it;
} }
else if (den > 0)
upperBound = std::min(upperBound, num/den);
else
lowerBound = std::max(lowerBound, num/den);
if (min < minLowBound) { ++it;
_minA(coeff, min);
return;
} }
min = min > 1 ? std::sqrt(min) : 1; if (lowerBound > upperBound)
//max = std::sqrt(max); _minAfast(coeff, mina);
else
mina = upperBound > 1 ? std::sqrt(upperBound) : 1;
} }
void MetricBasis::_minA(const fullMatrix<double> &coeff, double &mina) const void MetricBasis::_minAfast(const fullMatrix<double> &coeff, double &mina) const
{ {
double minq = coeff(0, 0); double minq = coeff(0, 0);
for (int i = 1; i < coeff.size1(); ++i) { for (int i = 1; i < coeff.size1(); ++i)
if (minq > coeff(i, 0)) minq = coeff(i, 0); minq = std::min(minq, coeff(i, 0));
}
double maxp = 0; double maxp = 0;
for (int i = 0; i < coeff.size1(); ++i) { for (int i = 0; i < coeff.size1(); ++i) {
...@@ -1311,65 +1539,146 @@ void MetricBasis::_minA(const fullMatrix<double> &coeff, double &mina) const ...@@ -1311,65 +1539,146 @@ void MetricBasis::_minA(const fullMatrix<double> &coeff, double &mina) const
maxp = std::max(maxp, tmp); maxp = std::max(maxp, tmp);
} }
mina = minq/maxp; mina = std::max(1., minq/maxp); // accept +inf as answer
if (mina < 1) mina = 1;
} }
void MetricBasis::_minK(const fullMatrix<double> &coeff, void MetricBasis::_minKfast(const fullMatrix<double> &coeff,
const fullVector<double> &jac, double &min) const const fullVector<double> &jac, double &minK) const
{ {
fullVector<double> r(coeff.size1()); fullVector<double> P(coeff.size1());
for (int i = 0; i < coeff.size1(); ++i) { for (int i = 0; i < coeff.size1(); ++i) {
r(i) = 0; P(i) = 0;
for (int l = 1; l < 7; ++l) { for (int l = 1; l < 7; ++l) {
r(i) += coeff(i, l) * coeff(i, l); P(i) += coeff(i, l) * coeff(i, l);
} }
r(i) = std::sqrt(r(i)); P(i) = std::sqrt(P(i));
} }
min = 1e10; double upperBound = std::numeric_limits<double>::infinity();
double lowerBound = -upperBound;
std::map<int, std::vector<IneqData> >::const_iterator itJ, itP; std::map<int, std::vector<IneqData> >::const_iterator itJ, itP;
itJ = _ineqJ2.begin(); itJ = _ineqJ2.begin();
itP = _ineqP3.begin(); itP = _ineqP3.begin();
if (_ineqP3.size() != _ineqJ2.size()) Msg::Fatal("sizes P3 %d, J2 %d", _ineqP3.size(), _ineqJ2.size()); // _ineqJ2 and _ineqP3 have the same size
//Msg::Warning("sizes %d %d", _ineqJ2.size(), _ineqP3.size()); while (itJ != _ineqJ2.end()) {
int count = 0; //if (itJ->first != itP->first) Msg::Fatal("not same hash %d %d", itJ->first, itP->first);
while (itJ != _ineqJ2.end() && itP != _ineqP3.end()) {
if (count >= (int)_ineqJ2.size()) Msg::Fatal("aaargh");
if (itJ->first != itP->first) Msg::Fatal("not same hash %d %d", itJ->first, itP->first);
double num = 0; double num = 0;
//Msg::Info("sizej %d", itJ->second.size()); for (unsigned int k = 0; k < itJ->second.size(); ++k) {
for (unsigned int l = 0; l < itJ->second.size(); ++l) { const int i = itJ->second[k].i;
const int i = itJ->second[l].i; const int j = itJ->second[k].j;
const int j = itJ->second[l].j; num += itJ->second[k].val * jac(i) * jac(j);
num += itJ->second[l].val * jac(i) * jac(j);
} }
double den = 0; double den = 0;
//Msg::Info("sizep %d", itP->second.size());
for (unsigned int l = 0; l < itP->second.size(); ++l) { for (unsigned int l = 0; l < itP->second.size(); ++l) {
const int i = itP->second[l].i; const int i = itP->second[l].i;
const int j = itP->second[l].j; const int j = itP->second[l].j;
const int k = itP->second[l].k; const int k = itP->second[l].k;
//Msg::Info("i%d j%d k%d", i, j, k); den += itP->second[l].val * P(i) * P(j) * P(k);
if (l>=itP->second.size()) Msg::Error("l %d/%d", l, itP->second.size());
if (i>=r.size() || j>=r.size()||k>=r.size() ) Msg::Fatal("i%d j%d k%d /%d (%dl%d)", i, j, k, r.size(), count, l);
den += itP->second[l].val * r(i) * r(j) * r(k);
} }
//Msg::Info("%g/%g = %g", num, den, num/den);
min = std::min(min, num/den); // minK has to satisfy: minK * den <= num
if (den == 0) {
if (num >= 0) {
++itJ;
++itP;
continue;
}
// otherwise, it's impossible to satisfy.
minK = 0;//TODO c'est mieux de retourner _minKfast(coeff, jac, minK); ?
return;
}
else if (den > 0)
upperBound = std::min(upperBound, num/den);
else
lowerBound = std::max(lowerBound, num/den);
++itJ; ++itJ;
++itP; ++itP;
++count;
} }
min = std::max(min, 0.);
if (lowerBound > upperBound)
minK = 0;//TODO c'est mieux de retourner _minKfast(coeff, jac, minK); ?
else
minK = std::max(.0, upperBound);
}
void MetricBasis::_minKsharp(const fullMatrix<double> &coeff,
const fullVector<double> &jac, double &minK) 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));
// fill only the half
for (int j = i; 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 upperBound = std::numeric_limits<double>::infinity();
double lowerBound = -upperBound;
std::map<int, std::vector<IneqData> >::const_iterator itJ, itP;
itJ = _ineqJ2.begin();
itP = _ineqP3.begin();
// _ineqJ2 and _ineqP3 have the same size
while (itJ != _ineqJ2.end()) {
double num = 0;
for (unsigned int k = 0; k < itJ->second.size(); ++k) {
const int i = itJ->second[k].i;
const int j = itJ->second[k].j;
num += itJ->second[k].val * jac(i) * jac(j);
}
double den = 0;
for (unsigned 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;
den += itP->second[l].val * 1/3*(Q(i,j)*P(k)+Q(i,k)*P(j)+Q(j,k)*P(i));
}
// minK has to satisfy: minK * den <= num
if (den == 0) {
if (num >= 0) {
++itJ;
++itP;
continue;
}
// otherwise, it's impossible to satisfy.
minK = 0;//TODO c'est mieux de retourner _minKfast(coeff, jac, minK); ?
return;
}
else if (den > 0)
upperBound = std::min(upperBound, num/den);
else
lowerBound = std::max(lowerBound, num/den);
++itJ;
++itP;
}
if (lowerBound > upperBound)
minK = 0;//TODO c'est mieux de retourner _minKfast(coeff, jac, minK); ?
else
minK = std::max(.0, upperBound);
} }
void MetricBasis::_maxAstKpos(const fullMatrix<double> &coeff, void MetricBasis::_maxAstKpos(const fullMatrix<double> &coeff,
const fullVector<double> &jac, double minK, double beta, double &maxa) const const fullVector<double> &jac, double minK, double beta, double &maxa) const
{ {
Msg::Fatal("Verifie si cette fonction est ok");
fullVector<double> P(coeff.size1()); fullVector<double> P(coeff.size1());
for (int i = 0; i < coeff.size1(); ++i) { for (int i = 0; i < coeff.size1(); ++i) {
P(i) = 0; P(i) = 0;
...@@ -1411,6 +1720,7 @@ void MetricBasis::_maxAstKpos(const fullMatrix<double> &coeff, ...@@ -1411,6 +1720,7 @@ void MetricBasis::_maxAstKpos(const fullMatrix<double> &coeff,
void MetricBasis::_maxAstKneg(const fullMatrix<double> &coeff, void MetricBasis::_maxAstKneg(const fullMatrix<double> &coeff,
const fullVector<double> &jac, double minK, double beta, double &maxa) const const fullVector<double> &jac, double minK, double beta, double &maxa) const
{ {
Msg::Fatal("Verifie si cette fonction est ok");
fullVector<double> P(coeff.size1()); fullVector<double> P(coeff.size1());
fullMatrix<double> Q(coeff.size1(), coeff.size1()); fullMatrix<double> Q(coeff.size1(), coeff.size1());
for (int i = 0; i < coeff.size1(); ++i) { for (int i = 0; i < coeff.size1(); ++i) {
...@@ -1462,6 +1772,7 @@ void MetricBasis::_maxAstKneg(const fullMatrix<double> &coeff, ...@@ -1462,6 +1772,7 @@ void MetricBasis::_maxAstKneg(const fullMatrix<double> &coeff,
void MetricBasis::_maxKstAfast(const fullMatrix<double> &coeff, void MetricBasis::_maxKstAfast(const fullMatrix<double> &coeff,
const fullVector<double> &jac, double mina, double beta, double &maxK) const const fullVector<double> &jac, double mina, double beta, double &maxK) const
{ {
Msg::Fatal("Verifie si cette fonction est ok");
fullVector<double> r(coeff.size1()); fullVector<double> r(coeff.size1());
for (int i = 0; i < coeff.size1(); ++i) { for (int i = 0; i < coeff.size1(); ++i) {
r(i) = 0; r(i) = 0;
...@@ -1503,6 +1814,7 @@ void MetricBasis::_maxKstAfast(const fullMatrix<double> &coeff, ...@@ -1503,6 +1814,7 @@ void MetricBasis::_maxKstAfast(const fullMatrix<double> &coeff,
void MetricBasis::_maxKstAsharp(const fullMatrix<double> &coeff, void MetricBasis::_maxKstAsharp(const fullMatrix<double> &coeff,
const fullVector<double> &jac, double mina, double beta, double &maxK) const const fullVector<double> &jac, double mina, double beta, double &maxK) const
{ {
Msg::Fatal("Verifie si cette fonction est ok");
fullVector<double> P(coeff.size1()); fullVector<double> P(coeff.size1());
fullMatrix<double> Q(coeff.size1(), coeff.size1()); fullMatrix<double> Q(coeff.size1(), coeff.size1());
for (int i = 0; i < coeff.size1(); ++i) { for (int i = 0; i < coeff.size1(); ++i) {
...@@ -1511,7 +1823,8 @@ void MetricBasis::_maxKstAsharp(const fullMatrix<double> &coeff, ...@@ -1511,7 +1823,8 @@ void MetricBasis::_maxKstAsharp(const fullMatrix<double> &coeff,
P(i) += coeff(i, l) * coeff(i, l); P(i) += coeff(i, l) * coeff(i, l);
} }
P(i) = std::sqrt(P(i)); P(i) = std::sqrt(P(i));
for (int j = 0; j < coeff.size1(); ++j) { // fill only the half
for (int j = i; j < coeff.size1(); ++j) {
Q(i, j) = 0; Q(i, j) = 0;
for (int l = 1; l < 7; ++l) { for (int l = 1; l < 7; ++l) {
Q(i, j) += coeff(i, l) * coeff(j, l); Q(i, j) += coeff(i, l) * coeff(j, l);
...@@ -1538,10 +1851,10 @@ void MetricBasis::_maxKstAsharp(const fullMatrix<double> &coeff, ...@@ -1538,10 +1851,10 @@ void MetricBasis::_maxKstAsharp(const fullMatrix<double> &coeff,
const int j = itP->second[l].j; const int j = itP->second[l].j;
const int k = itP->second[l].k; const int k = itP->second[l].k;
num += itP->second[l].val * coeff(i, 0) * coeff(j, 0) * coeff(k, 0); num += itP->second[l].val * coeff(i, 0) * coeff(j, 0) * coeff(k, 0);
if (j == k) if (j == k) // TODO verif ok
den += itP->second[l].val * Q(i,i) * P(i); den += itP->second[l].val * Q(i,i) * P(i);
else else
den += itP->second[l].val * 1/3*(Q(i,j)*P(k)+Q(i,k)*P(j)+Q(k,j)*P(i)); den += itP->second[l].val * 1/3*(Q(i,j)*P(k)+Q(i,k)*P(j)+Q(j,k)*P(i));
} }
min = std::min(min, num/den); min = std::min(min, num/den);
++itJ; ++itJ;
...@@ -1551,88 +1864,564 @@ void MetricBasis::_maxKstAsharp(const fullMatrix<double> &coeff, ...@@ -1551,88 +1864,564 @@ void MetricBasis::_maxKstAsharp(const fullMatrix<double> &coeff,
maxK = 1/beta*(mina*mina*mina-min); maxK = 1/beta*(mina*mina*mina-min);
} }
void MetricBasis::_computeTermBeta(double &a, double &K, void MetricBasis::_computeBoundBeta(const fullMatrix<double> &coeff,
double &dRda, double &term1, const fullVector<double> &jac,
double &phip) const double &beta, bool compLowBound) const
{
// compute a lower/upper bound for function beta = K/a^3
double upperBound = std::numeric_limits<double>::infinity();
double lowerBound = -upperBound;
// value in case of a failure
if (compLowBound) beta = 0;
else beta = 1;
std::map<int, std::vector<IneqData> >::const_iterator itJ, itP;
itJ = _ineqJ2.begin();
itP = _ineqP3.begin();
// _ineqJ2 and _ineqP3 have the same size
while (itJ != _ineqJ2.end()) {
double num = 0;
for (unsigned int k = 0; k < itJ->second.size(); ++k) {
const int i = itJ->second[k].i;
const int j = itJ->second[k].j;
num += itJ->second[k].val * jac(i) * jac(j);
}
double den = 0;
for (unsigned 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;
den += itP->second[l].val * coeff(i, 0) * coeff(j, 0) * coeff(k, 0);
}
// beta has to satisfy beta * den <= num if compLowBound = true
// or beta * den >= num if compLowBound = false
if (den == 0) {
if (compLowBound) {
if (num >= 0) {
++itJ;
++itP;
continue;
}
else return;
//TODO c'est mieux de retourner min(J^2)/max(q^3) ??
}
else {
if (num <= 0) {
++itJ;
++itP;
continue;
}
else return;
//TODO c'est mieux de retourner max(J^2)/min(q^3) ??
}
}
if (den > 0) {
if (compLowBound)
upperBound = std::min(upperBound, num/den);
else
lowerBound = std::max(lowerBound, num/den);
}
else {
if (compLowBound)
lowerBound = std::max(lowerBound, num/den);
else
upperBound = std::min(upperBound, num/den);
}
++itJ;
++itP;
}
if (lowerBound <= upperBound) {
if (compLowBound)
beta = std::max(.0, upperBound);
else
beta = std::min(1., lowerBound);
}
//TODO sinon, c'est mieux de retourner m..(J^2)/m..(q^3) ??
}
void MetricBasis::_computeBoundingCurve(const fullMatrix<double> &coeff,
const fullVector<double> &jac,
double &beta, double c,
bool compLowBound) const
{
// compute a lower/upper bounding curve K = beta * a^3 + c * a
// with c fixed.
//Msg::Info("BEGINNING %g %g", std::numeric_limits<double>::min(), std::numeric_limits<double>::denorm_min());
fullMatrix<double> Q(coeff.size1(), coeff.size1()); // (half filled!)
for (int i = 0; i < coeff.size1(); ++i) {
for (int j = i; 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 upperBound = std::numeric_limits<double>::infinity();
double lowerBound = -upperBound;
// values in case of a failure
if (compLowBound) beta = upperBound;
else beta = lowerBound;
std::map<int, std::vector<IneqData> >::const_iterator itJ, itP;
itJ = _ineqJ2.begin();
itP = _ineqP3.begin();
// _ineqJ2 and _ineqP3 have the same size
int kk = 0;
while (itJ != _ineqJ2.end()) {// && kk++ < 2) {// FIXME making tests
long double J2 = 0;
for (unsigned int k = 0; k < itJ->second.size(); ++k) {
const int i = itJ->second[k].i;
const int j = itJ->second[k].j;
J2 += itJ->second[k].val * jac(i) * jac(j);
//Msg::Info("(%d, %d) = (%g | %g, %g) -> J2:%Lg", i, j, itJ->second[k].val, jac(i), jac(j), J2);
}
long double q3 = 0, qp2 = 0;
for (unsigned 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;
q3 += itP->second[l].val * coeff(i, 0) * coeff(j, 0) * coeff(k, 0);
qp2 += itP->second[l].val * 1/3*(coeff(i, 0) * Q(j, k) +
coeff(j, 0) * Q(i, k) +
coeff(k, 0) * Q(i, j));
//Msg::Info("%d %d %d => %Lg %Lg", i, j, k, q3, qp2);
//Msg::Info("(%d %d %d) -> q3:%Lg qp2:%Lg", i, j, k, q3, qp2);
}
//Msg::Info("a%Lg K%Lg %Lg %Lg", std::sqrt(q3/qp2), std::sqrt(J2*J2*q3/qp2/qp2/qp2), std::sqrt(q3/qp2), std::sqrt(J2*J2*q3/qp2/qp2/qp2));
const long double num = J2 - ((long double)c)*qp2;
const long double den = q3;
// beta has to satisfy beta * den <= num if compLowBound = true
// or beta * den >= num if compLowBound = false
if (den == 0) {
if (compLowBound) {
if (num >= 0) {} // ok, nothing to do
else return;
//TODO c'est mieux de retourner min(..)/max(..) ??
}
else {
if (num <= 0) {} // ok, nothing to do
else return;
//TODO c'est mieux de retourner max(..)/min(..) ??
}
}
else if (den > 0) {
if (compLowBound)
upperBound = std::min(upperBound, (double)(num/den));
else
lowerBound = std::max(lowerBound, (double)(num/den));
}
else {
if (compLowBound)
lowerBound = std::max(lowerBound, (double)(num/den));
else
upperBound = std::min(upperBound, (double)(num/den));
}
//Msg::Info("in [%.15g, %.15g]", lowerBound, upperBound);
++itJ;
++itP;
}
if (lowerBound <= upperBound) {
if (compLowBound)
beta = upperBound;
else
beta = lowerBound;
}
//Msg::Info("ENDING");
//TODO sinon, c'est mieux de retourner m..(..)/m..(..) ??
}
void MetricBasis::_computeLowerBoundC(const fullMatrix<double> &coeff,
const fullVector<double> &jac,
double beta, double &c) const
{ {
double x0 = .5 * (K - a*a*a + 3*a); // compute a lower/upper bound of function C = K-beta*a^3
double sin, sqrt;
if (x0 > 1) { double upperBound = std::numeric_limits<double>::infinity();
const double p = -3; double lowerBound = -upperBound;
double q = -K + 2;
a = cubicCardanoRoot(p, q); // value in case of a failure
c = 0;
x0 = 1;
phip = M_PI / 3; fullVector<double> P(coeff.size1());
term1 = 1 + .5 * a; for (int i = 0; i < coeff.size1(); ++i) {
sin = std::sqrt(3.) / 2; P(i) = 0;
sqrt = 0; for (int l = 1; l < 7; ++l) {
} P(i) += coeff(i, l) * coeff(i, l);
else if (x0 < -1) {
K = -2 + a*a*a - 3*a;
x0 = -1;
phip = 2 * M_PI / 3;
term1 = 1 - .5 * a;
sin = std::sqrt(3.) / 2;
sqrt = 0;
} }
P(i) = std::sqrt(P(i));
}
std::map<int, std::vector<IneqData> >::const_iterator itJ, itP;
itJ = _ineqJ2.begin();
itP = _ineqP3.begin();
// _ineqJ2 and _ineqP3 have the same size
while (itJ != _ineqJ2.end()) {
double num = 0, den = 0;
for (unsigned 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 * P(i) * P(j) * P(k);
}
num = num * beta;
for (unsigned 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);
}
// c has to satisfy c * den <= num
if (den == 0) {
if (num >= 0) {
++itJ;
++itP;
continue;
}
else return;
//TODO c'est mieux de retourner min(J^2-beta q^3)/max(p^3) ?
}
if (den > 0)
upperBound = std::min(upperBound, num/den);
else { else {
phip = (std::acos(x0) + M_PI) / 3; lowerBound = std::max(lowerBound, num/den);
term1 = 1 + a * std::cos(phip); }
sin = std::sin(phip);
sqrt = std::sqrt(1-x0*x0); ++itJ;
++itP;
}
if (lowerBound <= upperBound) {
c = upperBound;
} }
dRda = sin * sqrt + .5 * term1 * (1-a*a); //TODO sinon, c'est mieux de retourner min(J^2-beta q^3)/max(p^3) ?
}
void MetricBasis::_computePseudoDerivatives(double a, double K,
double &dRda, double &dRdK)
{
// Return derivative without positive (but non-constant) coefficients
// Useful when only the sign of dRda and dRdK or the ratio dRda/dRdK
// are needed.
double w = _wSafe(a, K);
if (std::abs(w) > 1) {
Msg::Error("Cannot compute pseudo derivatives of R "
"outside the domain (w = %g)", w);
return;
}
const double phip = (std::acos(w) + M_PI) / 3;
const double C = 1 + a * std::cos(phip);
dRdK = C / 3;
dRda = 2 * std::sqrt(1-w*w) * std::sin(phip) + (1-a*a) * C;
}
void MetricBasis::_computeDerivatives(double a, double K,
double &dRda, double &dRdK,
double &dRdaa, double &dRdaK, double &dRdKK)
{
const double w = _wSafe(a, K);
if (std::abs(w) > 1) {
Msg::Error("Cannot compute derivatives of R outside the domain");
return;
}
const double phi = (std::acos(w)) / 3;
const double phip = phi + M_PI / 3;
const double S = 1./std::sqrt(1-w*w);
const double A = (1-a*a)*S;
const double sin = std::sin(phi);
const double sinp = std::sin(phip);
const double cosp = std::cos(phip);
const double C = 1 + a*cosp;
const double D = 1./(a+2*std::cos(phi));
static const double sq3 = std::sqrt(3);
const double coeff = sq3*D*D/18;
dRdK = coeff*6 * (C*S);
dRda = coeff*18 * (2*sinp + A*C);
dRdKK = coeff*S*S * (a*sinp - 4*C*sin*D + 3*w*C*S);
dRdaK = coeff*S*3 * (a*A*sinp + 3*w*A*C*S + 2*cosp - 4*C*(1+A*sin)*D);
dRdaa = coeff*9 * (3*w*A*A*C*S - 4*a*C*S + a*A*A*sinp - 4*(1+A*sin)*D*(2*sinp + A*C));
}
bool MetricBasis::_moveInsideDomain(double &a, double &K, bool bottomleftCorner)
{
// Note: For N-R (when moving 'a'), we compute the root of
// f(a) = .5 * (K - a^3 + 3*a) - 'w' where 'w' is -1 or 1.
double w = _w(a, K);
if (w > 1.) {
if (bottomleftCorner) {
// make sure to compute the right root (a>1) and
// avoid slopes that are too small:
a = std::max(a, 1.1);
const double tolerance = _tol/100;
// Note: N-R approaches the root from the right, we
while (w > 1 || w < 1-tolerance) {
double df = w - 1;
double slope = 1.5*(1-a*a);
a -= df/slope;
w = _w(a, K);
}
}
else {
K -= 2*w - 2;
}
return true;
}
else if (w < -1.) {
if (bottomleftCorner) {
K -= 2*w + 2;
}
else {
a = std::max(a, 2.);
const double tolerance = _tol/100;
while (std::abs(w+1) > tolerance) {
double df = w + 1;
double slope = 1.5*(1-a*a);
a -= df/slope;
w = _w(a, K);
}
}
return true;
}
return false;
}
double MetricBasis::_computeMinRAlongCurve(double beta, double c,
double mina, double maxa)
{
// Compute the minimum of R on the curve K = beta * a^3 + c
// where beta can be zero.
// 'mina' and 'maxa' stand for the limits of the curve. There must be at
// least 1 limit. To indicate that 'mina' or 'maxa' is not a limit, it is
// set to a negative value.
bool towardsPositivea;
double a, K, R;
if (mina < 0) {
a = maxa;
K = (beta*a*a + c)*a;
R = _R3Dsafe(a, K);
towardsPositivea = false;
}
else if (maxa < 0) {
a = mina;
K = (beta*a*a + c)*a;
R = _R3Dsafe(a, K);
towardsPositivea = true;
}
else {
// choose the more promising limit as starting point
double K1 = beta*mina*mina*mina + c;
double K2 = beta*maxa*maxa*maxa + c;
double R1 = _R3Dsafe(mina, K1);
double R2 = _R3Dsafe(maxa, K2);
if (R1 < R2) {
a = mina;
K = K1;
R = R1;
towardsPositivea = true;
}
else {
a = maxa;
K = K2;
R = R2;
towardsPositivea = false;
}
}
double da = 1, dK, dKK;
double dRda, dRdK, dRdaa, dRdaK, dRdKK;
double dRdC, dRdCC;
_computeDerivatives(a, K, dRda, dRdK, dRdaa, dRdaK, dRdKK);
dK = 3*beta*a*a+c;
dKK = 6*beta*a;
dRdC = dRda*da + dRdK*dK;
dRdCC = dRdaa*da*da + 2*dRdaK*da*dK + dRdKK*dK*dK + dRdK*dKK;
if ((towardsPositivea && dRdC > 0) || (!towardsPositivea && dRdC < 0)) {
Msg::Warning("Already the minimum");
return R;
}
while (.25*dRdC*dRdC/dRdCC > .01*_toleranceOnR(R)) {
a = a - dRdC / dRdCC;
K = (beta*a*a + c)*a;
R = _R3Dsafe(a, K);
_computeDerivatives(a, K, dRda, dRdK, dRdaa, dRdaK, dRdKK);
dK = 3*beta*a*a+c;
dKK = 6*beta*a;
dRdC = dRda*da + dRdK*dK;
dRdCC = dRdaa*da*da + 2*dRdaK*da*dK + dRdKK*dK*dK + dRdK*dKK;
}
/*Msg::Warning("Checking Newton Raphson solution... (%g, %g)", a, K);
double min = 1;
for (double aa = std::max(1., .5*a); aa < 2*a; aa += .001*a) {
double KK = (beta*aa*aa + c)*aa;
double w = _wSafe(aa, KK);
if (std::abs(w) <= 1) min = std::min(min, _R3Dsafe(aa, KK));
}
if (std::abs(min-R) < .01*_toleranceOnR(R))
Msg::Info("... aaand, it is good!");
else
Msg::Error("... and, it's a fail! %g vs %g (tol %g => %g)", R, min, _toleranceOnR(R), R-.01*_toleranceOnR(R));*/
return R;
}
bool MetricBasis::_intersectionCurveLeftCorner(double beta, double c,
double &a, double &K)
{
// on input, a and K are the bottom left corner
// on output, a and K are the intersection
// return true if the intersection is on the vertical
// We assume the intersection exists
// (i.e. the bounding curve is sufficiently sharp)
if (3*beta*a*a + c < 0) {
Msg::Error("The slope is negative, cannot compute the intersection");
return -1;
}
const double minK = K;
const double mina = a; //TODO remove (not needed)
K = (beta*a*a + c)*a;
if (K >= minK) return true;
// Find the root by Newton-Raphson
double dK = K-minK;
static double precision = std::numeric_limits<double>::epsilon() * 1e3;
while (std::abs(dK) > minK * precision) { //TODO meilleur approx
const double slope = 3*beta*a*a + c;
a -= dK / slope;
dK = (beta*a*a + c)*a-minK;
}
K = minK;
//Msg::Info("in %d it (%g, %g) -> (%g, %g) => diff(%g, %g)", k, mina, minK, a, K, a-mina, K-minK);
return false;
}
double MetricBasis::_computeAbscissaMinR(double aStart, double K)
{
double dRda, dRdK, dRdaa, dRdaK, dRdKK;
_computeDerivatives(aStart, K, dRda, dRdK, dRdaa, dRdaK, dRdKK);
double a = aStart;
while (dRda > MetricBasis::_tol*1e-3) {
double da = - dRda / dRdaa;
a += da;
_computeDerivatives(a, K, dRda, dRdK, dRdaa, dRdaK, dRdKK);
}
return a;
} }
double MetricBasis::_R3Dsafe(double q, double p, double J) double MetricBasis::_R3Dsafe(double q, double p, double J)
{ {
if (q > 1e5*p) { if (q >= p && p >= 0 && J == J &&
const double m = p*std::sqrt(3.)/q; p < std::numeric_limits<double>::infinity()) {
return (1-m) / (1+m); if (q == 0) return 0;
if (p == 0) return 1;
if (q == std::numeric_limits<double>::infinity()) {
Msg::Warning("q had infinite value in computation of 3D metric");
return 1;
} }
const double a = q/p; const double a = q/p;
const double K = J*J/p/p/p; const double K = J*J/p/p/p;
return _R3Dsafe(a, K); return _R3Dsafe(a, K);
} }
else {
Msg::Error("Wrong argument for 3d metric (%g, %g, %g)", q, p, J);
return -1;
}
}
double MetricBasis::_R3Dsafe(double a, double K) double MetricBasis::_R3Dsafe(double a, double K)
{ {
const double x = .5 * (K + (3 - a*a)*a); if (a >= 1 && K >= 0) {
if (x > 1+1e-7 || x < -1-1e-7) { // TODO speedup in function of _tol ?
Msg::Warning("x = %g (a,K) = (%g,%g)", x, a, K); // if (2/a < _tol) return 1.;
} // if (8/K < _tol*_tol*_tol) return 1.;
double ans; if (a == std::numeric_limits<double>::infinity())
if (x >= 1) ans = (a - 1) / (a + 2); return 1;
else if (x <= -1) ans = (a - 2) / (a + 1);
else { const double w = _wSafe(a, K);
const double phi = std::acos(x) / 3; if (std::abs(w) > 1-_tol/10) {
ans = (a + 2*std::cos(phi + 2*M_PI/3)) / (a + 2*std::cos(phi)); if (w > 0) return (a - 1) / (a + 2);
else return (a - 2) / (a + 1);
} }
if (ans < 0 || ans > 1) { const double phi = std::acos(w) / 3;
if (ans < 0) return 0; const double R = (a + 2*std::cos(phi + 2*M_PI/3)) / (a + 2*std::cos(phi));
else return 1; return std::max(.0, std::min(1., R));
}
else {
Msg::Error("Wrong argument for 3d metric (%g, %g)", a, K);
return -1;
} }
return ans;
} }
double MetricBasis::_R2Dsafe(double q, double p) double MetricBasis::_R2Dsafe(double q, double p)
{ {
if (q < 0 || p < 0 || p > q) if (q >= p && p >= 0 && p < std::numeric_limits<double>::infinity()) {
Msg::Error("wrong argument for 2d metric (%g, %g)", q, p); if (q == 0) return 0;
if (q == std::numeric_limits<double>::infinity()) {
Msg::Warning("q had infinite value in computation of 2D metric");
return 1;
}
return (q-p) / (q+p); return (q-p) / (q+p);
} }
else {
Msg::Error("Wrong argument for 2d metric (%g, %g)", q, p);
return -1;
}
}
double MetricBasis::_R2Dsafe(double a) double MetricBasis::_R2Dsafe(double a)
{ {
if (a < 1 if (a >= 1) {
#if !defined(_MSC_VER) if (a == std::numeric_limits<double>::infinity())
|| !std::isfinite(a) return 1;
#endif
)
Msg::Error("wrong argument for 2d metric (%g)", a);
return (a - 1) / (a + 1); return (a - 1) / (a + 1);
} }
else {
Msg::Error("Wrong argument for 2d metric (%g)", a);
return -1;
}
}
double MetricBasis::_toleranceOnR(double R)
{
const double sqRmax = std::max(.0, std::sqrt(R) - _tol);
return R - sqRmax*sqRmax;
}
...@@ -6,11 +6,13 @@ ...@@ -6,11 +6,13 @@
#ifndef _METRIC_BASIS_H_ #ifndef _METRIC_BASIS_H_
#define _METRIC_BASIS_H_ #define _METRIC_BASIS_H_
#include "JacobianBasis.h"
#include "fullMatrix.h"
#include <fstream> #include <fstream>
#include <cmath> #include <cmath>
#include "fullMatrix.h"
class MElement; class MElement;
class bezierBasis;
class JacobianBasis;
class GradientBasis;
class MetricBasis { class MetricBasis {
friend class MetricCoefficient; friend class MetricCoefficient;
...@@ -58,29 +60,46 @@ private: ...@@ -58,29 +60,46 @@ private:
public: public:
MetricBasis(int elementTag); MetricBasis(int elementTag);
// Get & Set methods
static void setTol(double tol) {_tol = tol;} static void setTol(double tol) {_tol = tol;}
const JacobianBasis* getJacobianForMetric() const {return _jacobian;} const JacobianBasis* getJacobianForMetric() const {return _jacobian;}
const bezierBasis* getBezier() const {return _bezier;} const bezierBasis* getBezier() const {return _bezier;}
int getNumMetricNodes() const {return _gradients->getNumSamplingPoints();} //int getNumMetricNodes() const {return _gradients->getNumSamplingPoints();}
// Compute min(R) with given tolerance
static double getMinR(MElement *el);
double computeMinR(MElement*) const;
// Sample R and return min
static double getMinSampledR(MElement *el, int order);
double computeMinSampledR(MElement*, int order) const;
// Compute R on corners and return min
static double getMinRCorner(MElement *el);
// TEST : sample min(r_min) and max(r_max) and return ratio of the two
static double getMinSampledGlobalRatio(MElement *el, int order);
double computeMinSampledGlobalRatio(MElement*, int order) const;
static double boundMinR(MElement *el); // TEST : sample R and return max
static double minSampledR(MElement *el, int order); static double getMaxSampledR(MElement *el, int order);
double getBoundMinR(MElement*) const; double computeMaxSampledR(MElement*, int order) const;
double getMinSampledR(MElement*, int order) const;
static double minRCorner(MElement *el);
//
template<bool ideal> template<bool ideal>
void getMetricCoeff(const fullMatrix<double> &nodes, void getMetricCoeff(const fullMatrix<double> &nodes,
fullMatrix<double> &coeff) const { fullMatrix<double> &coeff) const {
_fillCoeff<ideal>(_dim, _gradients, nodes, coeff); _fillCoeff<ideal>(_dim, _gradients, nodes, coeff);
} }
// void lag2Bez(const fullMatrix<double> &metCoeffLag,
// fullMatrix<double> &metCoeffBez) const {
// _bezier->matrixLag2Bez.mult(metCoeffLag, metCoeffBez);
// }
void lag2Bez(const fullMatrix<double> &metCoeffLag, void lag2Bez(const fullMatrix<double> &metCoeffLag,
fullMatrix<double> &metCoeffBez) const { fullMatrix<double> &metCoeffBez) const;
_bezier->matrixLag2Bez.mult(metCoeffLag, metCoeffBez);
}
// Metric order
static int metricOrder(int tag); static int metricOrder(int tag);
public: public:
...@@ -107,13 +126,13 @@ private: ...@@ -107,13 +126,13 @@ private:
void _fillInequalitiesPyr(int order); void _fillInequalitiesPyr(int order);
void _lightenInequalities(int&, int&, int&); //TODO change void _lightenInequalities(int&, int&, int&); //TODO change
void _computeRmin(const fullMatrix<double>&, const fullVector<double>&, void _computeBoundsRmin(const fullMatrix<double>&, const fullVector<double>&,
double &RminLag, double &RminBez) const; double &RminLag, double &RminBez) const;
void _computeRmax(const fullMatrix<double>&, const fullVector<double>&, void _computeBoundsRmax(const fullMatrix<double>&, const fullVector<double>&,
double &RmaxLag) const; double &RmaxLag) const;
void _getMetricData(const MElement*, MetricData*&) const; void _getMetricData(const MElement*, MetricData*&) const;
double _subdivideForRmin(MetricData*, double RminLag, double tol, MElement *el=NULL) const; double _subdivideForRmin(MetricData*, double &RminLag, double tol, MElement *el=NULL) const;
template<bool ideal> template<bool ideal>
static void _fillCoeff(int dim, const GradientBasis*, static void _fillCoeff(int dim, const GradientBasis*,
const fullMatrix<double> &nodes, fullMatrix<double> &coeff); const fullMatrix<double> &nodes, fullMatrix<double> &coeff);
...@@ -121,10 +140,9 @@ private: ...@@ -121,10 +140,9 @@ private:
const fullMatrix<double> &coeff, int num); const fullMatrix<double> &coeff, int num);
void _minMaxA(const fullMatrix<double>&, double &min, double &max) const; void _minMaxA(const fullMatrix<double>&, double &min, double &max) const;
void _minA(const fullMatrix<double>&, double &min) const; void _minAfast(const fullMatrix<double>&, double &min) const;
void _minK(const fullMatrix<double>&, const fullVector<double>&, double &min) const; void _minKfast(const fullMatrix<double>&, const fullVector<double>&, double &min) const;
void _computeTermBeta(double &a, double &K, double &dRda, void _minKsharp(const fullMatrix<double>&, const fullVector<double>&, double &min) const;
double &term1, double &phip) const;
void _maxAstKpos(const fullMatrix<double>&, const fullVector<double>&, void _maxAstKpos(const fullMatrix<double>&, const fullVector<double>&,
double minK, double beta, double &maxa) const; double minK, double beta, double &maxa) const;
void _maxAstKneg(const fullMatrix<double>&, const fullVector<double>&, void _maxAstKneg(const fullMatrix<double>&, const fullVector<double>&,
...@@ -133,11 +151,49 @@ private: ...@@ -133,11 +151,49 @@ private:
double mina, double beta, double &maxK) const; double mina, double beta, double &maxK) const;
void _maxKstAsharp(const fullMatrix<double>&, const fullVector<double>&, void _maxKstAsharp(const fullMatrix<double>&, const fullVector<double>&,
double mina, double beta, double &maxK) const; double mina, double beta, double &maxK) const;
void _computeBoundBeta(const fullMatrix<double>&,
const fullVector<double>&,
double &beta, bool lowerBound) const;
void _computeBoundingCurve(const fullMatrix<double>&,
const fullVector<double>&,
double &beta, double c, bool lowerBound) const;
void _computeLowerBoundC(const fullMatrix<double>&,
const fullVector<double>&,
double beta, double &c) const;
static bool _moveInsideDomain(double &a, double &K, bool bottomleftCorner);
static void _computePseudoDerivatives(double a, double K,
double &dRda, double &dRdK);
static void _computeDerivatives(double a, double K,
double &dRda, double &dRdK,
double &dRdaa, double &dRdaK, double &dRdKK);
static double _computeMinRAlongCurve(double beta, double c,
double mina, double maxa);
static bool _intersectionCurveLeftCorner(double beta, double c,
double &mina, double &minK);
static double _computeAbscissaMinR(double aStart, double K);
static double _R3Dsafe(double a, double K); static double _R3Dsafe(double a, double K);
static double _R3Dsafe(double q, double p, double J); static double _R3Dsafe(double q, double p, double J);
static double _R2Dsafe(double a); static double _R2Dsafe(double a);
static double _R2Dsafe(double q, double p); static double _R2Dsafe(double q, double p);
static inline double _w(double a, double K) {
return .5 * (K + (3-a*a)*a);
}
static inline double _wSafe(double a, double K) {
const double w = _w(a, K);
if (w > 1) {
if (w < 1+_tol/10) return 1;
else Msg::Error("outside the domain w(%g, %g) = %g", a, K, w);
}
else if (w < -1) {
if (w > -1-_tol/10) return -1;
else Msg::Error("outside the domain w(%g, %g) = %g", a, K, w);
}
return w;
}
static double _toleranceOnR(double R);
private: private:
class gterIneq { class gterIneq {
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment