From ceb735850325eee8e471bd8dd5320ec1f3275ccd Mon Sep 17 00:00:00 2001
From: Amaury Johnan <amjohnen@gmail.com>
Date: Fri, 28 Oct 2016 09:41:37 +0000
Subject: [PATCH] cleanup old class MetricBasis
---
Numeric/MetricBasis.cpp | 2160 ---------------------------------------
Numeric/MetricBasis.h | 200 ----
2 files changed, 2360 deletions(-)
delete mode 100644 Numeric/MetricBasis.cpp
delete mode 100644 Numeric/MetricBasis.h
diff --git a/Numeric/MetricBasis.cpp b/Numeric/MetricBasis.cpp
deleted file mode 100644
index b13874765b..0000000000
--- a/Numeric/MetricBasis.cpp
+++ /dev/null
@@ -1,2160 +0,0 @@
-// Gmsh - Copyright (C) 1997-2016 C. Geuzaine, J.-F. Remacle
-//
-// See the LICENSE.txt file for license information. Please report all
-// bugs and problems to the public mailing list <gmsh@onelab.info>.
-
-#include "MElement.h"
-#include "AnalyseCurvedMesh.h"
-#include "MetricBasis.h"
-#include "bezierBasis.h"
-#include "JacobianBasis.h"
-#include "BasisFactory.h"
-#include "pointsGenerators.h"
-#include "MHexahedron.h"
-#include "OS.h"
-#include <queue>
-#include <sstream>
-#include <limits>
-#include <iomanip>
-#include <cmath>
-
-double MetricBasis::_tol = 1e-3;
-
-//TODO Renommer fonctions plus explicitement (minmaxA,...) et rendre statique certaines fonctions
-
-namespace {
- /*
- double cubicCardanoRoot(double p, double q)
- {
- // solve the equation t^3 + p*t + q = 0
- double A = q*q/4 + p*p*p/27;
- if (A > 0) {
- double sq = std::sqrt(A);
- 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 {
- double module = std::sqrt(-p*p*p/27);
- double ang = std::acos(-q/2/module);
- return 2 * std::pow(module, 1/3.) * std::cos(ang/3);
- }
- }
- */
- int nChoosek(int n, int k)
- {
- if (n < k || k < 0) {
- Msg::Error("Wrong argument for combination. n %d k %d", n, k);
- return 1;
- }
-
- if (k > n/2) k = n-k;
- if (k == 1)
- return n;
- if (k == 0)
- return 1;
-
- int c = 1;
- for (int i = 1; i <= k; i++, n--) (c *= n) /= i;
- return c;
- }
-
- double symRand(double f = 1)
- {
- return f * (rand()%2001 - 1000) / 1000.;
- }
-}
-
-MetricBasis::MetricBasis(int tag) :
- _jacobian(NULL), _type(ElementType::ParentTypeFromTag(tag)),
- _dim(ElementType::DimensionFromTag(tag))
-{
- //Msg::Fatal("Verifie si cette fonction est ok 1");
- const bool serendip = false;
- const int metOrder = metricOrder(tag);
- const int jacOrder = 3*metOrder/2;
-
- // get bezier and gradients for metric space
- FuncSpaceData *metricSpace = NULL;
- if (_type != TYPE_PYR)
- metricSpace = new FuncSpaceData(true, tag, metOrder, &serendip);
- else
- metricSpace = new FuncSpaceData(true, tag, false, metOrder+2,
- metOrder, &serendip);
- _gradients = BasisFactory::getGradientBasis(*metricSpace);
- _bezier = BasisFactory::getBezierBasis(*metricSpace);
- delete metricSpace;
-
- // get jacobian
- FuncSpaceData *jacSpace = NULL;
- if (_type == TYPE_TET || _type == TYPE_HEX || _type == TYPE_PRI) {
- jacSpace = new FuncSpaceData(true, tag, jacOrder, &serendip);
- }
- 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,
- jacOrder, &serendip);
- }
- else if (_type == TYPE_TRI || _type == TYPE_QUA) {
- // jacobian not needed
- }
- else
- Msg::Fatal("metric not implemented for element tag %d", tag);
-
- if (jacSpace) {
- _jacobian = BasisFactory::getJacobianBasis(*jacSpace);
- delete jacSpace;
- }
-
- _fillInequalities(metOrder);
-}
-
-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());
- double m = metric->computeMinR(el);
- //Msg::Info("returning %g", m);
- return m;
-}
-
-double MetricBasis::computeMinR(MElement *el) const
-{
- int nSampPnts = _gradients->getNumSamplingPoints();
- int nMapping = _gradients->getNumMapNodes();
-
- // Nodes
- fullMatrix<double> nodes(nMapping, 3);
- el->getNodesCoord(nodes);
-
- // Jacobian coefficients
- fullVector<double> *jac = NULL;
- if (_jacobian) {
- fullVector<double> jacLag(_jacobian->getNumJacNodes());
- jac = new fullVector<double>(_jacobian->getNumJacNodes());
- _jacobian->getSignedIdealJacobian(nodes, jacLag);
- _jacobian->lag2Bez(jacLag, *jac);
- }
-
- // Metric coefficients
- fullMatrix<double> metCoeffLag;
- _fillCoeff<true>(el->getDim(), _gradients, nodes, metCoeffLag);
- fullMatrix<double> *metCoeff;
- metCoeff = new fullMatrix<double>(nSampPnts, metCoeffLag.size2());
- _bezier->matrixLag2Bez.mult(metCoeffLag, *metCoeff);
-
- /*static double inf = std::numeric_limits<double>::infinity();
- 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::getMinSampledR(MElement *el, int order)
-{
- MetricBasis *metric = (MetricBasis*)BasisFactory::getMetricBasis(el->getTypeForMSH());
- return metric->computeMinSampledR(el, order);
-}
-
-double MetricBasis::computeMinSampledR(MElement *el, int deg) const
-{
- //Msg::Fatal("Verifie si cette fonction est ok b");
- 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 min = R(0, 1);
- for (int i = 1; i < R.size1(); ++i)
- min = std::min(min, R(i, 1));
-
- delete md;
- return min;
-}
-
-double MetricBasis::getMaxSampledR(MElement *el, int order)
-{
- MetricBasis *metric = (MetricBasis*)BasisFactory::getMetricBasis(el->getTypeForMSH());
- 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
- fullMatrix<double> nodes(nMapping, 3);
- el->getNodesCoord(nodes);
-
- // Jacobian coefficients
- fullVector<double> jacLag(jacobian->getNumJacNodes());
- jacobian->getSignedJacobian(nodes, jacLag);
-
- // Metric coefficients
- fullMatrix<double> metCoeffLag;
- _fillCoeff<false>(el->getDim(), gradients, nodes, metCoeffLag);
-
- // Compute min_corner(R)
- return _computeMinlagR(jacLag, metCoeffLag, nSampPnts);
-}
-
-void MetricBasis::_fillInequalities(int metricOrder)
-{
- //Msg::Fatal("Verifie si cette fonction est ok");
- if (_type == TYPE_PYR) {
- _fillInequalitiesPyr(metricOrder);
- return;
- }
-
- fullMatrix<int> exp(_bezier->_exponents.size1(), _bezier->_exponents.size2());
- for (int i = 0; i < _bezier->_exponents.size1(); ++i) {
- for (int j = 0; j < _bezier->_exponents.size2(); ++j) {
- exp(i, j) = static_cast<int>(_bezier->_exponents(i, j) + .5);
- }
- }
-
- int ncoeff = _gradients->getNumSamplingPoints();
- int dimSimplex = _bezier->getDimSimplex();
- int dim = _bezier->getDim();
-
- int countP3 = 0, countJ2 = 0, countA = 0;
- for (int i = 0; i < ncoeff; i++) {
- for (int j = i; j < ncoeff; j++) {
- double num = 1, den = 1;
- {
- int compl1 = metricOrder;
- int compl2 = metricOrder;
- int compltot = 2*metricOrder;
- for (int k = 0; k < dimSimplex; k++) {
- num *= nChoosek(compl1, exp(i, k))
- * nChoosek(compl2, exp(j, k));
- den *= nChoosek(compltot, exp(i, k) + exp(j, k));
- compl1 -= exp(i, k);
- compl2 -= exp(j, k);
- compltot -= exp(i, k) + exp(j, k);
- }
- for (int k = dimSimplex; k < dim; k++) {
- num *= nChoosek(metricOrder, exp(i, k))
- * nChoosek(metricOrder, exp(j, k));
- den *= nChoosek(2*metricOrder, exp(i, k) + exp(j, k));
- }
- }
-
- if (i != j) num *= 2;
-
- ++countA;
- int hash = 0;
- for (int l = 0; l < dim; l++) {
- hash += (exp(i, l)+exp(j, l)) * pow_int(2*metricOrder+1, l);
- }
- _ineqA[hash].push_back(IneqData(num/den, i, j));
-
- for (int k = j; k < ncoeff; ++k) {
- double num = 1, den = 1;
- {
- int compl1 = metricOrder;
- int compl2 = metricOrder;
- int compl3 = metricOrder;
- int compltot = 3*metricOrder;
- for (int l = 0; l < dimSimplex; l++) {
- num *= nChoosek(compl1, exp(i, l))
- * nChoosek(compl2, exp(j, l))
- * nChoosek(compl3, exp(k, l));
- den *= nChoosek(compltot, exp(i, l) + exp(j, l) + exp(k, l));
- compl1 -= exp(i, l);
- compl2 -= exp(j, l);
- compl3 -= exp(k, l);
- compltot -= exp(i, l) + exp(j, l) + exp(k, l);
- }
- for (int l = dimSimplex; l < dim; l++) {
- num *= nChoosek(metricOrder, exp(i, l))
- * nChoosek(metricOrder, exp(j, l))
- * nChoosek(metricOrder, exp(k, l));
- den *= nChoosek(3*metricOrder, exp(i, l) + exp(j, l) + exp(k, l));
- }
- }
-
- if (i == j) {
- if (j != k) num *= 3;
- }
- else {
- if (j == k || i == k) {
- num *= 3;
- }
- else num *= 6;
- }
-
- ++countP3;
- int hash = 0;
- 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 (_dim == 2) {
- _lightenInequalities(countJ2, countP3, countA);
- return;
- }
-
- exp.resize(_jacobian->getBezier()->_exponents.size1(),
- _jacobian->getBezier()->_exponents.size2());
- for (int i = 0; i < _jacobian->getBezier()->_exponents.size1(); ++i) {
- for (int j = 0; j < _jacobian->getBezier()->_exponents.size2(); ++j) {
- exp(i, j) = static_cast<int>(_jacobian->getBezier()->_exponents(i, j) + .5);
- }
- }
- int njac = _jacobian->getNumJacNodes();
- for (int i = 0; i < njac; i++) {
- for (int j = i; j < njac; j++) {
- int order = metricOrder/2*3;
- double num = 1, den = 1;
- {
- int compl1 = order;
- int compl2 = order;
- int compltot = 2*order;
- for (int k = 0; k < dimSimplex; k++) {
- num *= nChoosek(compl1, exp(i, k))
- * nChoosek(compl2, exp(j, k));
- den *= nChoosek(compltot, exp(i, k) + exp(j, k));
- compl1 -= exp(i, k);
- compl2 -= exp(j, k);
- compltot -= exp(i, k) + exp(j, k);
- }
- }
- for (int k = dimSimplex; k < dim; k++) {
- num *= nChoosek(order, exp(i, k))
- * nChoosek(order, exp(j, k));
- den *= nChoosek(2*order, exp(i, k) + exp(j, k));
- }
-
- if (i != j) num *= 2;
-
- ++countJ2;
- int hash = 0;
- for (int k = 0; k < dim; k++) {
- hash += (exp(i, k)+exp(j, k)) * pow_int(2*order+1, k);
- }
- _ineqJ2[hash].push_back(IneqData(num/den, i, j));
- }
- }
-
- _lightenInequalities(countJ2, countP3, countA);
-}
-
-void MetricBasis::_fillInequalitiesPyr(int metricOrder)
-{
- //Msg::Fatal("Verifie si cette fonction est ok");
- fullMatrix<int> exp(_bezier->_exponents.size1(), _bezier->_exponents.size2());
- for (int i = 0; i < _bezier->_exponents.size1(); ++i) {
- for (int j = 0; j < _bezier->_exponents.size2(); ++j) {
- exp(i, j) = static_cast<int>(_bezier->_exponents(i, j) + .5);
- }
- }
-
- int ncoeff = _gradients->getNumSamplingPoints();
-
- int countP3 = 0, countJ2 = 0, countA = 0;
- for (int i = 0; i < ncoeff; i++) {
- for (int j = i; j < ncoeff; j++) {
- double num, den;
- num = nChoosek(metricOrder+2, exp(i, 0))
- * nChoosek(metricOrder+2, exp(i, 1))
- * nChoosek(metricOrder , exp(i, 2))
- * nChoosek(metricOrder+2, exp(j, 0))
- * nChoosek(metricOrder+2, exp(j, 1))
- * nChoosek(metricOrder , exp(j, 2));
- den = nChoosek(2*metricOrder+4, exp(i, 0) + exp(j, 0))
- * nChoosek(2*metricOrder+4, exp(i, 1) + exp(j, 1))
- * nChoosek(2*metricOrder , exp(i, 2) + exp(j, 2));
-
- if (i != j) num *= 2;
-
- ++countA;
- int hash = 0;
- for (int l = 0; l < 3; l++) {
- hash += (exp(i, l)+exp(j, l)) * pow_int(2*metricOrder+1, l);
- }
- _ineqA[hash].push_back(IneqData(num/den, i, j));
-
-
- for (int k = j; k < ncoeff; ++k) {
- double num, den;
- num = nChoosek(metricOrder+2, exp(i, 0))
- * nChoosek(metricOrder+2, exp(i, 1))
- * nChoosek(metricOrder , exp(i, 2))
- * nChoosek(metricOrder+2, exp(j, 0))
- * nChoosek(metricOrder+2, exp(j, 1))
- * nChoosek(metricOrder , exp(j, 2))
- * nChoosek(metricOrder+2, exp(k, 0))
- * nChoosek(metricOrder+2, exp(k, 1))
- * nChoosek(metricOrder , exp(k, 2));
- den = nChoosek(3*metricOrder+6, exp(i, 0) + exp(j, 0) + exp(k, 0))
- * nChoosek(3*metricOrder+6, exp(i, 1) + exp(j, 1) + exp(k, 1))
- * nChoosek(3*metricOrder , exp(i, 2) + exp(j, 2) + exp(k, 2));
-
- if (i == j) {
- if (j != k) num *= 3;
- }
- else {
- if (j == k || i == k) {
- num *= 3;
- }
- else num *= 6;
- }
-
- ++countP3;
- int hash = 0;
- for (int l = 0; l < 3; 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));
- }
- }
- }
-
- exp.resize(_jacobian->getBezier()->_exponents.size1(),
- _jacobian->getBezier()->_exponents.size2());
- for (int i = 0; i < _jacobian->getBezier()->_exponents.size1(); ++i) {
- for (int j = 0; j < _jacobian->getBezier()->_exponents.size2(); ++j) {
- exp(i, j) = static_cast<int>(_jacobian->getBezier()->_exponents(i, j) + .5);
- }
- }
- int njac = _jacobian->getNumJacNodes();
- for (int i = 0; i < njac; i++) {
- for (int j = i; j < njac; j++) {
- int order = metricOrder/2*3;
-
- double num, den;
- num = nChoosek(order+3, exp(i, 0))
- * nChoosek(order+3, exp(i, 1))
- * nChoosek(order , exp(i, 2))
- * nChoosek(order+3, exp(j, 0))
- * nChoosek(order+3, exp(j, 1))
- * nChoosek(order , exp(j, 2));
- den = nChoosek(2*order+6, exp(i, 0) + exp(j, 0))
- * nChoosek(2*order+6, exp(i, 1) + exp(j, 1))
- * nChoosek(2*order , exp(i, 2) + exp(j, 2));
-
- if (i != j) num *= 2;
-
- ++countJ2;
- int hash = 0;
- for (int k = 0; k < 3; k++) {
- hash += (exp(i, k)+exp(j, k)) * pow_int(2*order+1, k);
- }
- _ineqJ2[hash].push_back(IneqData(num/den, i, j));
- }
- }
-
- _lightenInequalities(countJ2, countP3, countA);
-}
-
-void MetricBasis::_lightenInequalities(int &countj, int &countp, int &counta)
-{
- //Msg::Fatal("Verifie si cette fonction est ok");
- double tol = .0;
- std::map<int, std::vector<IneqData> >::iterator it, itbeg[3], itend[3];
-
- int cnt[3] = {0,0,0};
- itbeg[0] = _ineqJ2.begin();
- itbeg[1] = _ineqP3.begin();
- itbeg[2] = _ineqA.begin();
- itend[0] = _ineqJ2.end();
- itend[1] = _ineqP3.end();
- itend[2] = _ineqA.end();
- for (int k = 0; k < 3; ++k) {
- it = itbeg[k];
- while (it != itend[k]) {
- std::sort(it->second.begin(), it->second.end(), gterIneq());
-
- double rmved = .0;
- while (it->second.size() && rmved + it->second.back().val <= tol) {
- rmved += it->second.back().val;
- it->second.pop_back();
- ++cnt[k];
- }
- const double factor = 1-rmved;
- for (unsigned int i = 0; i < it->second.size(); ++i) {
- it->second[i].val /= factor;
- }
- ++it;
- }
- }
- countj -= cnt[0];
- countp -= cnt[1];
- counta -= cnt[2];
-}
-
-bool MetricBasis::validateBezierForMetricAndJacobian()
-{
- Msg::Fatal("Verifie si cette fonction est ok f");
- Msg::Info("Testing Bezier interpolation and subdivision "
- "for jacobien and metric on all element types...");
- 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;
- }
-#if defined(HAVE_PLUGINS)
- GMSH_AnalyseCurvedMeshPlugin plugin = GMSH_AnalyseCurvedMeshPlugin();
- plugin.test(elements, numElem, dim);
-#endif
- }
-
-
- //=================
- return true; //987;
- //=================
-
- /*Msg::Info("Testing Bezier interpolation and subdivision "
- "for jacobien and metric on all element types...");
- int numError = 0;
-
- // Parameters
- const int numType = 6;
- const int acceptedTypes[numType] = {TYPE_TRI, TYPE_QUA, TYPE_TET,
- TYPE_PYR, TYPE_PRI, TYPE_HEX};
- const int maxOrder = 3; // at least 3 (so that serendip tet are tested)
- const int numElem = 5; // at least 2 (first is reference element)
- const int numSubdiv = 2; // at least 1
- const int numSampPnt = 10; // at least 1
- const double toleranceTensor = 1e-11; // at most 1e-5 (metric takes values in [0,1])
- double tolerance; // computed in function of tag
-
- //
- static const double epsilon = std::numeric_limits<double>::epsilon();
-
- for (int tag = 1; tag <= MSH_NUM_TYPE; ++tag) {
- if (tag > 66 && tag < 71) continue; //not conventional elements
- if (tag > 75 && tag < 79) continue; //no element tag 76, 77, 78...
-
- // Check if accepted type
- const int type = ElementType::ParentTypeFromTag(tag);
- bool knownType = false;
- for (int i = 0; i < numType; ++i) {
- knownType = (knownType || type == acceptedTypes[i]);
- }
- if (!knownType) continue;
-
- const int order = ElementType::OrderFromTag(tag);
- const int dim = ElementType::DimensionFromTag(tag);
- const bool serendip = ElementType::SerendipityFromTag(tag) > 1;
-
- // Skip p0 elements and elements for which order > 'maxOrder'
- // and skip for now serendipty pyramids (not implemented)
- if (order < 1 || order > maxOrder) continue;
- if (type == TYPE_PYR && serendip) continue;
-
- Msg::Info("... testing element tag %d", tag);
-
- // Compute tolerance
- tolerance = epsilon * pow_int(10, order*dim);
- if (type == TYPE_PYR) tolerance = std::max(tolerance, epsilon*1e9);
-
- // 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;
- std::vector<MVertex*> vertices(nodes.size1());
- for (int i = 0; i < nodes.size1(); ++i) {
- vertices[i] = new MVertex(nodes(i, 0) + symRand(range),
- dim > 1 ? nodes(i, 1) + symRand(range) : 0,
- dim > 2 ? nodes(i, 2) + symRand(range) : 0);
- }
- MElement *el = MElement::createElement(tag, vertices);
- if (!el) {
- Msg::Error("MElement was unable to create element for tag %d", tag);
- ++numError;
- }
-
- numError += validateBezierForMetricAndJacobian(el, numSampPnt, numSubdiv,
- toleranceTensor, tolerance);
- }
- }
-
- if (numError) Msg::Error("Validation of Bezier terminated with %d errors!", numError);
- else Msg::Info("Validation of Bezier terminated without errors");
- return numError;*/
-}
-
-int MetricBasis::validateBezierForMetricAndJacobian(MElement *el,
- int numSampPnt,
- int numSubdiv,
- double toleranceTensor,
- double tolerance)
-{
- Msg::Fatal("Verifie si cette fonction est ok g");
- const int tag = el->getTypeForMSH();
- const MetricBasis *metricBasis = BasisFactory::getMetricBasis(tag);
-
- // compare the two metric
- fullMatrix<double> metric_Bez(numSampPnt, 2);
- fullVector<int> isub(numSubdiv);
- fullMatrix<double> uvw(numSampPnt, 3);
- metricBasis->interpolateAfterNSubdivisions(el, numSubdiv, numSampPnt,
- isub, uvw, metric_Bez);
-
- int numBadMatch = 0;
- int numBadMatchTensor = 0;
- double maxBadMatch = 0;
- double maxBadMatchTensor = 0;
- for (int isamp = 0; isamp < numSampPnt; ++isamp) {
- double dum[3];
- double &u = uvw(isamp, 0);
- double &v = uvw(isamp, 1);
- double &w = uvw(isamp, 2);
- double metric_Lag = el->getEigenvaluesMetric(u, v, w, dum);
-
- double diff = std::abs(metric_Lag - metric_Bez(isamp, 0));
- double diffTensor = std::abs(metric_Lag - metric_Bez(isamp, 1));
-
- if (diffTensor > toleranceTensor) {
- ++numBadMatchTensor;
- maxBadMatchTensor = std::max(maxBadMatchTensor, diffTensor);
- }
- else if (diff > tolerance) {
- ++numBadMatch;
- maxBadMatch = std::max(maxBadMatch, diff);
- }
- }
-
- if (numBadMatchTensor > .2*numSampPnt) {
- Msg::Error("Too much errors "
- "even when computing by metric tensor (max %g)", maxBadMatchTensor);
- return 1;
- }
- else if (numBadMatch > .5*numSampPnt) {
- Msg::Error("Too much errors (max %g)", maxBadMatch);
- return 1;
- }
- return 0;
-}
-
-void MetricBasis::interpolate(const MElement *el,
- const MetricData *md,
- const fullMatrix<double> &nodes,
- fullMatrix<double> &R) const
-{
- //Msg::Fatal("Verifie si cette fonction est ok h");
- fullMatrix<double> &metcoeffs = *md->_metcoeffs, *metric = new fullMatrix<double>;
- fullVector<double> &jaccoeffs = *md->_jaccoeffs, *jac = new fullVector<double>;
-
- _bezier->interpolate(metcoeffs, nodes, *metric);
-
- R.resize(nodes.size1(), 4);
-
- switch (metcoeffs.size2()) {
- case 1:
- for (int i = 0; i < R.size1(); ++i)
- R(i, 0) = R(i, 1) = 1;
- break;
-
- 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) {
- // Compute from q, p
- double p = pow((*metric)(i, 1), 2);
- p += pow((*metric)(i, 2), 2);
- p = std::sqrt(p);
- R(i, 0) = std::sqrt(_R2Dsafe((*metric)(i, 0), p));
- // Comppute from tensor
- fullMatrix<double> metricTensor(2, 2);
- metricTensor(0, 0) = (*metric)(i, 0) + (*metric)(i, 1);
- metricTensor(1, 1) = (*metric)(i, 0) - (*metric)(i, 1);
- metricTensor(0, 1) = metricTensor(1, 0) = (*metric)(i, 2);
- fullVector<double> valReal(2), valImag(2);
- fullMatrix<double> vecLeft(2, 2), vecRight(2, 2);
- metricTensor.eig(valReal, valImag, vecLeft, vecRight, true);
- 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;
-
- case 7:
- ((MetricBasis*)this)->file << "a K myR eigR" << std::endl;
- ((MetricBasis*)this)->file << _bezier->getNumLagCoeff() << std::endl;
- _jacobian->getBezier()->interpolate(jaccoeffs, nodes, *jac);
- for (int i = 0; i < R.size1(); ++i) {
- // Compute from q, p, J
- double p = 0;
- for (int k = 1; k < 7; ++k) p += pow((*metric)(i, k), 2);
- p = std::sqrt(p);
- R(i, 0) = std::sqrt(_R3Dsafe((*metric)(i, 0), p, (*jac)(i)));
- // Compute from tensor
- fullMatrix<double> metricTensor(3, 3);
- for (int k = 0; k < 3; ++k) {
- static double fact1 = std::sqrt(6.);
- static double fact2 = std::sqrt(3.);
- const int ki = k%2;
- const int kj = std::min(k+1, 2);
- metricTensor(k, k) = (*metric)(i, k+1)*fact1 + (*metric)(i, 0);
- metricTensor(ki, kj) = metricTensor(kj, ki) = (*metric)(i, k+4)*fact2;
- }
- fullVector<double> valReal(3), valImag(3);
- fullMatrix<double> vecLeft(3, 3), vecRight(3, 3);
- metricTensor.eig(valReal, valImag, vecLeft, vecRight, true);
- 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;
-
- default:
- Msg::Error("Wrong number of functions for metric: %d",
- metcoeffs.size2());
- }
-
- delete jac;
- delete metric;
-}
-
-void MetricBasis::interpolateAfterNSubdivisions(
- const MElement *el, int numSubdiv, int numPnt,
- fullVector<int> &isub,
- fullMatrix<double> &uvw,
- fullMatrix<double> &metric) const
-{
- Msg::Fatal("Verifie si cette fonction est ok i");
- // Interpolate metric after 'numSub' random subdivision at
- // 'numPnt' random points.
- // Return: isub, the subdomain tag of each subdivision,
- // uvw, the reference points at which metric has been interpolated.
- // metric, the interpolation.
-
- MetricData *md;
- _getMetricData(el, md);
-
- // Keep trace of subdomain to be able to compute uvw.
- // (Ensure to have the tag for the complete element):
- const nodalBasis *mapBasis = BasisFactory::getNodalBasis(el->getTypeForMSH());
- fullMatrix<double> nodes;
- mapBasis->getReferenceNodesForBezier(nodes);
-
- const int nSub = _bezier->getNumDivision();
- const int numCoeff = md->_metcoeffs->size2();
- const int numMetPnts = md->_metcoeffs->size1();
- const int numJacPnts = md->haveJac() ? md->_jaccoeffs->size() : 0;
- const int numNodPnts = nodes.size1();
-
- const bezierBasis *bezierMapping;
- if (el->getType() != TYPE_PYR)
- bezierMapping = BasisFactory::getBezierBasis(el->getTypeForMSH());
- else {
- FuncSpaceData data(true, el->getTypeForMSH(), false,
- el->getPolynomialOrder(), el->getPolynomialOrder());
- bezierMapping = BasisFactory::getBezierBasis(data);
- }
-
- for (int k = 0; k < numSubdiv; ++k) {
- fullMatrix<double> subcoeffs, subnodes;
- fullVector<double> subjac;
- _bezier->subdivideBezCoeff(*md->_metcoeffs, subcoeffs);
- bezierMapping->subdivideBezCoeff(nodes, subnodes);
-
- if (_dim == 3)
- _jacobian->getBezier()->subdivideBezCoeff(*md->_jaccoeffs, subjac);
- delete md;
-
- isub(k) = std::rand() % nSub;
-
- fullMatrix<double> *coeff = new fullMatrix<double>(numMetPnts, numCoeff);
- coeff->copy(subcoeffs, isub(k) * numMetPnts, numMetPnts, 0, numCoeff, 0, 0);
- nodes.copy(subnodes, isub(k) * numNodPnts, numNodPnts, 0, _dim, 0, 0);
- fullVector<double> *jac = NULL;
- if (_dim == 3) {
- jac = new fullVector<double>(numJacPnts);
- jac->copy(subjac, isub(k) * numJacPnts, numJacPnts, 0);
- }
-
- md = new MetricData(coeff, jac);
- }
-
- // compute a random convex combination of reference nodes
- fullMatrix<double> subuvw(uvw.size1(), _dim);
- {
- int tagPrimary = ElementType::getTag(el->getType(), 1);
- const nodalBasis *primMapBasis = BasisFactory::getNodalBasis(tagPrimary);
- fullMatrix<double> refNodes;
- primMapBasis->getReferenceNodes(refNodes);
- double *c = new double[refNodes.size1()];
-
- for (int k = 0; k < uvw.size1(); ++k) {
- double sum = 0;
- int exp = 1 + std::rand() % 5;
- for (int i = 0; i < refNodes.size1(); ++i) {
- c[i] = pow_int((std::rand() % 1000) / 1000., exp);
- sum += c[i];
- }
- for (int i = 0; i < refNodes.size1(); ++i) {
- c[i] /= sum;
- subuvw(k, 0) += c[i] * refNodes(i, 0);
- if (_dim > 1) subuvw(k, 1) += c[i] * refNodes(i, 1);
- if (_dim > 2) subuvw(k, 2) += c[i] * refNodes(i, 2);
- }
- }
-
- delete[] c;
- }
-
- interpolate(el, md, subuvw, metric);
- bezierMapping->interpolate(nodes, subuvw, uvw, false);
- delete md;
-}
-
-void MetricBasis::statsForMatlab(MElement *el, int deg, MetricData *md) const
-{
- //return;
- //Msg::Fatal("Verifie si cette fonction est ok h");
- fullMatrix<double> samplingPoints;
- bool serendip = false;
- gmshGeneratePoints(FuncSpaceData(el, deg, &serendip), samplingPoints);
-
- if (!md) _getMetricData(el, md);
-
- static unsigned int aa = 0;
- if (md->_num < 100000 && ++aa < 1000) {
- std::stringstream name;
- name << "metricStat_" << 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().c_str(), std::fstream::out);
-
- {
- ((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;
- fullVector<double> *jac = md->_jaccoeffs;
- _minMaxA(*coeff, mina, maxa);
-
- if (el->getDim() == 3) {
- _minKfast(*coeff, *jac, minKf);
- _minKsharp(*coeff, *jac, minKs);
- double minK = minKf;
- _moveInsideDomain(mina, minK, true);
- double p_dRda, p_dRdK;
- _computePseudoDerivatives(mina, minK, p_dRda, p_dRdK);
- double slope = -p_dRda/p_dRdK;
- c = .5*(3*minK/mina - slope);
- _computeBoundingCurve(*coeff, *jac, beta_m, c, true);
- _computeBoundingCurve(*coeff, *jac, beta_M, c, false);
- beta = -p_dRda/(p_dRdK*3*mina*mina);
- if (beta < 0)
- _computeLowerBoundC(*coeff, *jac, beta, c_m);
- else
- c_m = -1;
- }
- else {
- minKs = minKf = c = beta_m = beta_M = beta = c_m = -1;
- }
-
- ((MetricBasis*)this)->file << std::setprecision(15);
- ((MetricBasis*)this)->file << mina << " " << minKf << " " << minKs << " ";
- ((MetricBasis*)this)->file << c << " " << beta_m << " " << beta_M << " ";
- ((MetricBasis*)this)->file << beta << " " << c_m << std::endl;
- ((MetricBasis*)this)->file << std::setprecision(8);
- }
- }
-
- fullMatrix<double> R;
- interpolate(el, md, samplingPoints, R);
- ((MetricBasis*)this)->file.close();
-
- return;
-
- if (R.size1() < 1) {
- ((MetricBasis*)this)->file << -1 << " ";
- ((MetricBasis*)this)->file.close();
- return;
- }
-
- double min = R(0, 1);
- for (int i = 1; i < R.size1(); ++i)
- min = std::min(min, R(i, 1));
-
- //((MetricBasis*)this)->file << min << " ";
- ((MetricBasis*)this)->file.close();
- return;
-}
-
-int MetricBasis::metricOrder(int tag)
-{
- const int parentType = ElementType::ParentTypeFromTag(tag);
- const int order = ElementType::OrderFromTag(tag);
-
- switch (parentType) {
- case TYPE_PNT : return 0;
-
- case TYPE_LIN : return order;
-
- case TYPE_TRI :
- case TYPE_TET :
- case TYPE_PYR : return 2*order-2;
-
- case TYPE_QUA :
- case TYPE_PRI :
- case TYPE_HEX : return 2*order;
-
- default :
- Msg::Error("Unknown element type %d, return order 0", parentType);
- return 0;
- }
-}
-
-void MetricBasis::_computeBoundsRmin(
- const fullMatrix<double> &coeff, const fullVector<double> &jac,
- double &sqrtRminLag, double &sqrtRminBez) const
-{
- sqrtRminLag = _computeMinlagR(jac, coeff, _bezier->getNumLagCoeff());
- if (sqrtRminLag <= _tol) {
- sqrtRminBez = .0;
- return;
- }
-
- if (coeff.size2() == 3) { // 2d element
- double mina, dummy;
- _minMaxA(coeff, mina, dummy);
- sqrtRminBez = std::sqrt(_R2Dsafe(mina));
- return;
- }
-
- double minK, mina, dummy;
- _minKfast(coeff, jac, minK);
- _minMaxA(coeff, mina, dummy);
-
- _moveInsideDomain(mina, minK, true);
-
- double p_dRda, p_dRdK;
- _computePseudoDerivatives(mina, minK, p_dRda, p_dRdK);
-
- if (p_dRda < 0) {
- //Msg::Info("case 1 & 2");
- // 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);
- {
- double beta = (minK/mina-c)/mina/mina;
- if (std::abs(p_dRda + p_dRdK * (3*beta*mina*mina+c)) > 1e-12) {
- Msg::Error("Derivative along curve not zero %g", p_dRda + p_dRdK * (3*beta*mina*mina+c));
- }
- }
-
- double a_int = mina, K_int = minK;
- bool onVertical = _intersectionCurveLeftCorner(beta, c, a_int, K_int);
- if (onVertical) {
- // then the minimum is on the curve, we prefer to return this
- sqrtRminBez = std::sqrt(_R3Dsafe(mina, minK));
- //Msg::Info("onvertical");
- return;
- }
- /*if (a_int - mina < ???) {
- sqrtRminBez = std::sqrt(_R3Dsafe(a_int, minK));
- return;
- }*/
-
- bool haveToCompute_a0;
- if (_moveInsideDomain(a_int, K_int, true)) {
- haveToCompute_a0 = true;
- }
- else {
- double p_dRda, p_dRdK;
- _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);
- }
- haveToCompute_a0 = p_dRda > 0;
- }
-
- if (haveToCompute_a0) {
- double a0 = _computeAbscissaMinR(mina, minK);
- sqrtRminBez = std::sqrt(_R3Dsafe(a0, minK));
- }
- else {
- sqrtRminBez = std::sqrt(_R3Dsafe(a_int, K_int));
- }
- return;
- }
- else if (p_dRdK < 0) {
- //Msg::Info("case 4 & 5");
- double slope = -p_dRda/p_dRdK;
- double c = .5*(3*minK/mina - slope);
- double beta;
- _computeBoundingCurve(coeff, jac, beta, c, false);
- double a_int = mina, K_int = minK;
-
- if (_intersectionCurveLeftCorner(beta, c, a_int, K_int)) {
- /*if (std::abs(a_int - mina) < _tol / 100 &&
- 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 {
- if (_moveInsideDomain(a_int, K_int, false))
- sqrtRminBez = std::sqrt(_R3Dsafe(a_int, K_int));
- else
- sqrtRminBez = std::sqrt(_computeMinRAlongCurve(beta, c, a_int, -1));
- }
- return;
- }
- else {
- //Msg::Info("case 3");
- sqrtRminBez = std::sqrt(_R3Dsafe(mina, minK));
- return;
- }
-}
-
-void MetricBasis::_computeBoundsRmax(
- const fullMatrix<double> &coeff, const fullVector<double> &jac,
- double &RmaxLag) const
-{
- Msg::Fatal("Verifie si cette fonction est ok k");
- RmaxLag = 0.;
-
- 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);
- const double a = q/p;
- if (a > 1e4) {
- RmaxLag = std::max(RmaxLag, std::sqrt((a - std::sqrt(3.)) / (a + std::sqrt(3.))));
- }
- else {
- const double tmpR = _R3Dsafe(a, jac(i)/p/p*jac(i)/p);
- RmaxLag = std::max(RmaxLag, std::sqrt(tmpR));
- }
- }
-}
-
-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;
- const bool for3d = md->_jaccoeffs;
- const int numCoeff = md->_metcoeffs->size2();
- const int numMetPnts = md->_metcoeffs->size1();
- const int numJacPnts = for3d ? md->_jaccoeffs->size() : 0;
- const int numSub = _bezier->getNumDivision();
- subdomains.push(md);
-
- std::vector<fullVector<double>*> trash;
-
- int k = 0;
- while (RminLag - subdomains.top()->_RminBez > tol && subdomains.size() < 25000) {
- fullMatrix<double> *subcoeffs, *coeff;
- fullVector<double> *subjac, *jac = NULL;
-
- MetricData *current = subdomains.top();
- subcoeffs = new fullMatrix<double>(numSub*numMetPnts, numCoeff);
- _bezier->subDivisor.mult(*current->_metcoeffs, *subcoeffs);
- if (for3d) {
- subjac = new fullVector<double>(numSub*numJacPnts);
- _jacobian->getBezier()->subDivisor.mult(*current->_jaccoeffs, *subjac);
- }
- int depth = current->_depth;
- int num = current->_num;
- delete current;
- subdomains.pop();
-
- for (int i = 0; i < numSub; ++i) {
- coeff = new fullMatrix<double>(numMetPnts, numCoeff);
- coeff->copy(*subcoeffs, i * numMetPnts, numMetPnts, 0, numCoeff, 0, 0);
- if (for3d) {
- jac = new fullVector<double>;
- jac->setAsProxy(*subjac, i * numJacPnts, numJacPnts);
- }
- double minLag, minBez;
- _computeBoundsRmin(*coeff, *jac, minLag, minBez);
- RminLag = std::min(RminLag, minLag);
- int newNum = num + (i+1) * pow_int(10, depth);
- if (depth > 8) newNum = 999999999;
- MetricData *metData = new MetricData(coeff, jac, minBez, depth+1, newNum);
- //statsForMatlab(el, 15, metData);
- //if (el) statsForMatlab(el, 20, metData);
- subdomains.push(metData);
- }
- if (for3d) trash.push_back(subjac);
- delete subcoeffs;
- ++k;
- }
-
- //Msg::Info("num subdiv %d", k);
-
- double ans = subdomains.top()->_RminBez;
- while (subdomains.size()) {
- md = subdomains.top();
- subdomains.pop();
- delete md;
- }
- for (unsigned int i = 0; i < trash.size(); ++i) {
- delete trash[i];
- }
-
- return ans;
-}
-
-void MetricBasis::_getMetricData(const MElement *el, MetricData *&md) const
-{
- int nSampPnts = _gradients->getNumSamplingPoints();
- int nMapping = _gradients->getNumMapNodes();
-
- // Nodes
- fullMatrix<double> nodes(nMapping, 3);
- el->getNodesCoord(nodes);
-
- // Jacobian coefficients
- fullVector<double> *jac = NULL;
- if (_dim == 3) {
- fullVector<double> jacLag(_jacobian->getNumJacNodes());
- jac = new fullVector<double>(_jacobian->getNumJacNodes());
- _jacobian->getSignedIdealJacobian(nodes, jacLag);
- _jacobian->lag2Bez(jacLag, *jac);
- }
-
- // Metric coefficients
- fullMatrix<double> metCoeffLag;
- _fillCoeff<true>(el->getDim(), _gradients, nodes, metCoeffLag);
- fullMatrix<double> *metCoeff;
- metCoeff = new fullMatrix<double>(nSampPnts, metCoeffLag.size2());
- _bezier->matrixLag2Bez.mult(metCoeffLag, *metCoeff);
-
- md = new MetricData(metCoeff, jac, -1, 0, 0);
-}
-
-template<bool ideal>
-void MetricBasis::_fillCoeff(int dim, const GradientBasis *gradients,
- const fullMatrix<double> &nodes, fullMatrix<double> &coeff)
-{
- //Msg::Fatal("Verifie si cette fonction est ok o");
- const int nSampPnts = gradients->getNumSamplingPoints();
-
- switch (dim) {
- case 0 :
- return;
- case 1 :
- Msg::Fatal("Should not be here, metric for 1d not implemented");
- break;
-
- case 2 :
- {
- fullMatrix<double> dxydX(nSampPnts,3), dxydY(nSampPnts,3);
- if (ideal)
- gradients->getIdealGradientsFromNodes(nodes, &dxydX, &dxydY, NULL);
- else
- gradients->getGradientsFromNodes(nodes, &dxydX, &dxydY, NULL);
-
- coeff.resize(nSampPnts, 3);
- for (int i = 0; i < nSampPnts; i++) {
- const double &dxdX = dxydX(i,0), &dydX = dxydX(i,1);
- const double &dxdY = dxydY(i,0), &dydY = dxydY(i,1);
- double dzdX = 0, dzdY = 0;
- if (nodes.size2() > 2) {
- dzdX = dxydX(i,2);
- dzdY = dxydY(i,2);
- }
- const double dvxdX = dxdX*dxdX + dydX*dydX + dzdX*dzdX;
- const double dvxdY = dxdY*dxdY + dydY*dydY + dzdY*dzdY;
- coeff(i, 0) = (dvxdX + dvxdY) / 2;
- coeff(i, 1) = dvxdX - coeff(i, 0);
- coeff(i, 2) = (dxdX*dxdY + dydX*dydY + dzdX*dzdY);
- }
- }
- break;
-
- case 3 :
- {
- fullMatrix<double> dxyzdX(nSampPnts,3), dxyzdY(nSampPnts,3), dxyzdZ(nSampPnts,3);
- if (ideal)
- gradients->getIdealGradientsFromNodes(nodes, &dxyzdX, &dxyzdY, &dxyzdZ);
- else
- gradients->getGradientsFromNodes(nodes, &dxyzdX, &dxyzdY, &dxyzdZ);
-
- coeff.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;
- coeff(i, 0) = (dvxdX + dvxdY + dvxdZ) / 3;
- static double fact1 = 1./std::sqrt(6.);
- static double fact2 = 1./std::sqrt(3.);
- coeff(i, 1) = fact1 * (dvxdX - coeff(i, 0));
- coeff(i, 2) = fact1 * (dvxdY - coeff(i, 0));
- coeff(i, 3) = fact1 * (dvxdZ - coeff(i, 0));
- coeff(i, 4) = fact2 * (dxdX*dxdY + dydX*dydY + dzdX*dzdY);
- coeff(i, 5) = fact2 * (dxdZ*dxdY + dydZ*dydY + dzdZ*dzdY);
- coeff(i, 6) = fact2 * (dxdX*dxdZ + dydX*dydZ + dzdX*dzdZ);
- }
- }
- break;
- }
-}
-
-double MetricBasis::_computeMinlagR(const fullVector<double> &jac,
- const fullMatrix<double> &coeff, int num)
-{
- double Rmin = 1.;
-
- switch (coeff.size2()) {
- case 7:
- for (int i = 0; i < num; ++i) {
- if (jac(i) <= 0.) return 0;
-
- const 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);
-
- Rmin = std::min(Rmin, _R3Dsafe(q, p, jac(i)));
- }
- break;
-
- case 3:
- for (int i = 0; i < num; ++i) {
- const double &q = coeff(i, 0);
- double p = pow_int(coeff(i, 1), 2) + pow_int(coeff(i, 2), 2);
- p = std::sqrt(p);
-
- Rmin = std::min(Rmin, _R2Dsafe(q, p));
- }
- break;
-
- default:
- Msg::Error("coeff have not right number of column");
- return -1;
- }
-
- return std::sqrt(Rmin);
-}
-
-void MetricBasis::_minMaxA(const fullMatrix<double> &coeff,
- double &mina, double &max) const
-{
- // TODO calculer le max ?
- double upperBound = std::numeric_limits<double>::infinity();
- double lowerBound = -upperBound;
- std::map<int, std::vector<IneqData> >::const_iterator it = _ineqA.begin();
-
- while (it != _ineqA.end()) {
- double num = 0;
- double den = 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 < coeff.size2(); ++l) {
- tmp += coeff(i, l) * coeff(j, l);
- }
- den += it->second[k].val * tmp;
- num += it->second[k].val * coeff(i, 0) * coeff(j, 0);
- }
-
- // mina has to satisfy: mina^2 * den <= num
- if (den == 0) {
- if (num >= 0) {} // ok, nothing to do
- else {// impossible to satisfy.
- _minAfast(coeff, mina);
- return;
- }
- }
- else if (den > 0)
- upperBound = std::min(upperBound, num/den);
- else
- lowerBound = std::max(lowerBound, num/den);
-
- ++it;
- }
-
- if (lowerBound > upperBound)
- _minAfast(coeff, mina);
- else
- mina = upperBound > 1 ? std::sqrt(upperBound) : 1;
-}
-
-void MetricBasis::_minAfast(const fullMatrix<double> &coeff, double &mina) const
-{
- double minq = coeff(0, 0);
- for (int i = 1; i < coeff.size1(); ++i)
- minq = std::min(minq, coeff(i, 0));
-
- double maxp = 0;
- for (int i = 0; i < coeff.size1(); ++i) {
- double tmp = 0;
- for (int j = 1; j < 7; ++j) {
- tmp += pow_int(coeff(i, j), 2);
- }
- maxp = std::max(maxp, tmp);
- }
-
- mina = std::max(1., minq/maxp); // accept +inf as answer
-}
-
-void MetricBasis::_minKfast(const fullMatrix<double> &coeff,
- const fullVector<double> &jac, double &minK) const
-{
- fullVector<double> P(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));
- }
-
- 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()) {
- //if (itJ->first != itP->first) Msg::Fatal("not same hash %d %d", itJ->first, itP->first);
- 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 * P(i) * P(j) * P(k);
- }
-
- // 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::_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::_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
-{
- // compute a lower/upper bound of function C = K-beta*a^3
-
- double upperBound = std::numeric_limits<double>::infinity();
- double lowerBound = -upperBound;
-
- // value in case of a failure
- c = 0;
-
- fullVector<double> P(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));
- }
-
- 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 {
- lowerBound = std::max(lowerBound, num/den);
- }
-
- ++itJ;
- ++itP;
- }
-
- if (lowerBound <= upperBound) {
- c = upperBound;
- }
- //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 false;
- }
-
- 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;
- int n = 0;
- while (std::abs(dK) > minK * precision && ++n < 100) { //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)
-{
- if (q >= p && p >= 0 && J == J &&
- p < std::numeric_limits<double>::infinity()) {
- 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 K = J*J/p/p/p;
- 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)
-{
- if (a >= 1 && K >= 0) {
- // TODO speedup in function of _tol ?
- // if (2/a < _tol) return 1.;
- // if (8/K < _tol*_tol*_tol) return 1.;
-
- if (a == std::numeric_limits<double>::infinity())
- return 1;
-
- const double w = _wSafe(a, K);
- if (std::abs(w) > 1-_tol/10) {
- if (w > 0) return (a - 1) / (a + 2);
- else return (a - 2) / (a + 1);
- }
-
- const double phi = std::acos(w) / 3;
- const double R = (a + 2*std::cos(phi + 2*M_PI/3)) / (a + 2*std::cos(phi));
- return std::max(.0, std::min(1., R));
- }
- else {
- Msg::Error("Wrong argument for 3d metric (%g, %g)", a, K);
- return -1;
- }
-}
-
-double MetricBasis::_R2Dsafe(double q, double p)
-{
- if (q >= p && p >= 0 && p < std::numeric_limits<double>::infinity()) {
- 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);
- }
- else {
- Msg::Error("Wrong argument for 2d metric (%g, %g)", q, p);
- return -1;
- }
-}
-
-double MetricBasis::_R2Dsafe(double a)
-{
- if (a >= 1) {
- if (a == std::numeric_limits<double>::infinity())
- return 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;
-}
diff --git a/Numeric/MetricBasis.h b/Numeric/MetricBasis.h
deleted file mode 100644
index 5d0a389f96..0000000000
--- a/Numeric/MetricBasis.h
+++ /dev/null
@@ -1,200 +0,0 @@
-// Gmsh - Copyright (C) 1997-2016 C. Geuzaine, J.-F. Remacle
-//
-// See the LICENSE.txt file for license information. Please report all
-// bugs and problems to the public mailing list <gmsh@onelab.info>.
-
-#ifndef _METRIC_BASIS_H_
-#define _METRIC_BASIS_H_
-
-#include <fstream>
-#include <cmath>
-#include "fullMatrix.h"
-class MElement;
-class bezierBasis;
-class JacobianBasis;
-class GradientBasis;
-
-class MetricBasis {
- friend class MetricCoefficient;
- friend class GMSH_AnalyseCurvedMeshPlugin;
-
-private:
- const JacobianBasis *_jacobian;
- const GradientBasis *_gradients;
- const bezierBasis *_bezier;
- static double _tol;
- const int _type, _dim;
-
- 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) {}
- };
-
- std::map<int, std::vector<IneqData> > _ineqJ2, _ineqP3, _ineqA;
-
- class MetricData {
- public:
- fullMatrix<double> *_metcoeffs;
- fullVector<double> *_jaccoeffs;
- double _RminBez;
- int _depth, _num;
-
- public:
- MetricData(fullMatrix<double> *m, fullVector<double> *j,
- double r = -1, int d = -1, int num = -1) :
- _metcoeffs(m), _jaccoeffs(j), _RminBez(r), _depth(d), _num(num) {}
- ~MetricData() {
- delete _metcoeffs;
- delete _jaccoeffs;
- }
-
- bool haveJac() {return _jaccoeffs != NULL;}
- };
-
-public:
- MetricBasis(int elementTag);
-
- // Get & Set methods
- static void setTol(double tol) {_tol = tol;}
- const JacobianBasis* getJacobianForMetric() const {return _jacobian;}
- const bezierBasis* getBezier() const {return _bezier;}
- //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;
-
- // TEST : sample R and return max
- static double getMaxSampledR(MElement *el, int order);
- double computeMaxSampledR(MElement*, int order) const;
-
-
- //
- template<bool ideal>
- void getMetricCoeff(const fullMatrix<double> &nodes,
- fullMatrix<double> &coeff) const {
- _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,
- fullMatrix<double> &metCoeffBez) const;
-
- // Metric order
- static int metricOrder(int tag);
-
-public:
- // Validation for computation of Bezier coefficients & subdivision
- // of Jacobian determinant and Metric stuffs
- static bool validateBezierForMetricAndJacobian();
- static int validateBezierForMetricAndJacobian(MElement *el,
- int numSampPnt,
- int numSubdiv,
- double toleranceTensor,
- double tolerance);
- void statsForMatlab(MElement *el, int deg, MetricData *md) const;
- void interpolate(const MElement*, const MetricData*,
- const fullMatrix<double> &nodes, fullMatrix<double> &R) const;
- void interpolateAfterNSubdivisions(const MElement *el,
- int numSubdiv, int numPnt,
- fullVector<int> &isub,
- fullMatrix<double> &uvw,
- fullMatrix<double> &metric) const;
-
-private:
- void _fillInequalities(int order);
- void _fillInequalitiesPyr(int order);
- void _lightenInequalities(int&, int&, int&); //TODO change
-
- void _computeBoundsRmin(const fullMatrix<double>&, const fullVector<double>&,
- double &RminLag, double &RminBez) const;
- void _computeBoundsRmax(const fullMatrix<double>&, const fullVector<double>&,
- double &RmaxLag) const;
- void _getMetricData(const MElement*, MetricData*&) const;
-
- double _subdivideForRmin(MetricData*, double &RminLag, double tol, MElement *el=NULL) const;
- template<bool ideal>
- static void _fillCoeff(int dim, const GradientBasis*,
- const fullMatrix<double> &nodes, fullMatrix<double> &coeff);
- static double _computeMinlagR(const fullVector<double> &jac,
- const fullMatrix<double> &coeff, int num);
-
- void _minMaxA(const fullMatrix<double>&, double &min, double &max) const;
- void _minAfast(const fullMatrix<double>&, double &min) const;
- void _minKfast(const fullMatrix<double>&, const fullVector<double>&, double &min) const;
- void _minKsharp(const fullMatrix<double>&, const fullVector<double>&, double &min) 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 q, double p, double J);
- static double _R2Dsafe(double a);
- 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:
- 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;
- }
- };
-};
-
-#endif
--
GitLab