diff --git a/Numeric/JacobianBasis.cpp b/Numeric/JacobianBasis.cpp
index 24157004ccf46ff9f94923bf06db392c62557027..62a4a0e778cb1d04f6ad9b125a55b02875c3a532 100644
--- a/Numeric/JacobianBasis.cpp
+++ b/Numeric/JacobianBasis.cpp
@@ -88,11 +88,15 @@ void GradientBasis::getGradientsFromNodes(const fullMatrix<double> &nodes,
if (dxyzdZ) gradShapeMatZ.mult(nodes, *dxyzdZ);
}
-JacobianBasis::JacobianBasis(int tag)
+JacobianBasis::JacobianBasis(int tag, int jacOrder)
{
const int parentType = ElementType::ParentTypeFromTag(tag);
const int order = ElementType::OrderFromTag(tag);
- const int jacobianOrder = JacobianBasis::jacobianOrder(parentType, order);
+ int jacobianOrder;
+ if (jacOrder < 0)
+ jacobianOrder = JacobianBasis::jacobianOrder(parentType, order);
+ else
+ jacobianOrder = jacOrder;
const int primJacobianOrder = JacobianBasis::jacobianOrder(parentType, 1);
fullMatrix<double> lagPoints; // Sampling points
diff --git a/Numeric/JacobianBasis.h b/Numeric/JacobianBasis.h
index cdcd1adab34bc8bed90e26ccce24618af7a0d65c..aec12be7e6996f1bec5fe4d9825afa0b674ab7c9 100644
--- a/Numeric/JacobianBasis.h
+++ b/Numeric/JacobianBasis.h
@@ -44,7 +44,7 @@ class JacobianBasis {
public :
const bezierBasis *bezier;
- JacobianBasis(int tag);
+ JacobianBasis(int tag, int jacOrder = -1);
// Get methods
inline int getNumJacNodes() const { return numJacNodes; }
diff --git a/Numeric/MetricBasis.cpp b/Numeric/MetricBasis.cpp
index 1b96d6f6df1a7e6f6bee81011092d8dad4375a57..1f3b277870b08d14c10326c8071db0dfc1bcaa09 100644
--- a/Numeric/MetricBasis.cpp
+++ b/Numeric/MetricBasis.cpp
@@ -8,24 +8,56 @@
#include "BasisFactory.h"
#include "pointsGenerators.h"
#include <cmath>
-
+#include <queue>
+#include "OS.h"
#if defined(HAVE_SCIP)
-#include <scip/scip.h>
-#include <soplex.h>
+#include "scip/scip.h"
+#include "scip/scipdefplugins.h"
+#include "soplex/soplex.h"
+#include <sstream>
+#include "scip/cons.h"
+#include "scip/struct_scip.h"
+#endif
-void MaFonctionScip()
-{
- SCIP *scip;
- SCIPcreate(&scip);
+namespace {
+ double cubicCardanoRoot(double p, double q)
+ {
+ double sq = std::sqrt(q*q/4 + p*p*p/27);
+ return std::pow(-q/2+sq, 1/3.) + std::pow(-q/2-sq, 1/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;
+ }
}
-#endif
MetricCoefficient::MetricCoefficient(MElement *el) : _element(el)
{
const int tag = el->getTypeForMSH();
const int type = ElementType::ParentTypeFromTag(tag);
const int metricOrder = MetricCoefficient::metricOrder(tag);
- _jacobian = BasisFactory::getJacobianBasis(tag);
+ if (type == TYPE_HEX || type == TYPE_PRI) {
+ int order = ElementType::OrderFromTag(tag);
+ _jacobian = new JacobianBasis(tag, 3*order);
+ }
+ else if (type == TYPE_TET)
+ _jacobian = BasisFactory::getJacobianBasis(tag);
+ else
+ Msg::Fatal("metric not implemented for element tag %d", tag);
_gradients = BasisFactory::getGradientBasis(tag, metricOrder);
_bezier = BasisFactory::getBezierBasis(type, metricOrder);
@@ -67,13 +99,400 @@ MetricCoefficient::MetricCoefficient(MElement *el) : _element(el)
_coefficientsLag(i, 1) = dvxdX - _coefficientsLag(i, 0);
_coefficientsLag(i, 2) = dvxdY - _coefficientsLag(i, 0);
_coefficientsLag(i, 3) = dvxdZ - _coefficientsLag(i, 0);
- _coefficientsLag(i, 4) = dxdX*dxdY + dydX*dydY + dzdX*dzdY;
- _coefficientsLag(i, 5) = dxdZ*dxdY + dydZ*dydY + dzdZ*dzdY;
- _coefficientsLag(i, 6) = dxdX*dxdZ + dydX*dydZ + dzdX*dzdZ;
+ const double fact = std::sqrt(2);
+ _coefficientsLag(i, 4) = fact * (dxdX*dxdY + dydX*dydY + dzdX*dzdY);
+ _coefficientsLag(i, 5) = fact * (dxdZ*dxdY + dydZ*dydY + dzdZ*dzdY);
+ _coefficientsLag(i, 6) = fact * (dxdX*dxdZ + dydX*dydZ + dzdX*dzdZ);
}
}
break;
}
+
+ _fillInequalities(metricOrder);
+}
+
+void MetricCoefficient::_fillInequalities(int metricOrder)
+{
+ std::map<int, int> result2vect;
+ std::vector<std::pair<int,int> > v;
+ int dimSimplex = _bezier->_dimSimplex;
+ int dim = _bezier->getDim();
+ fullMatrix<double> exp = _bezier->_exponents;
+ int ncoeff = _coefficientsLag.size1();
+ Msg::Info("ncoeff %d", ncoeff);
+ int countP3 = 0, countJ2 = 0;
+ double tol = .99;
+ double tol2 = .001;
+ for (int i = 0; i < ncoeff; i++) {
+ for (int j = i; j < ncoeff; j++) {
+ double A = 1, B = 1;
+ {
+ int compl1 = metricOrder;
+ int compl2 = metricOrder;
+ int compl3 = 2*metricOrder;
+ for (int k = 0; k < dimSimplex; k++) {
+ A *= nChoosek(compl1, (int) exp(i, k))
+ * nChoosek(compl2, (int) exp(j, k));
+ B *= nChoosek(compl3, (int) exp(i, k) + (int) exp(j, k));
+ compl1 -= (int) exp(i, k);
+ compl2 -= (int) exp(j, k);
+ compl3 -= (int) exp(i, k) + (int) exp(j, k);
+ }
+ }
+ for (int k = dimSimplex; k < dim; k++) {
+ A *= nChoosek(metricOrder, (int) exp(i, k))
+ * nChoosek(metricOrder, (int) exp(j, k));
+ B *= nChoosek(2*metricOrder, (int) exp(i, k) + (int) exp(j, k));
+ }
+
+ if (i != j) A *= 2;
+
+ if (A > tol*B) {
+ v.push_back(std::make_pair(i, j));
+ int ind;
+ int hash = 0;
+ for (int k = 0; k < dim; k++) {
+ hash += (int) (exp(i, k)+exp(j, k)) * pow_int(2*metricOrder+1, k);
+ }
+ std::map<int, int>::iterator it;
+ if ((it = result2vect.find(hash)) != result2vect.end()) {
+ ind = it->second;
+ }
+ else {
+ ind = _ineq2.size();
+ _ineq2.push_back(std::map<std::pair<int, int>, double>());
+ result2vect[hash] = ind;
+ }
+ _ineq2[ind][std::make_pair(i, j)] = A/B;
+ }
+ for (int k = j; k < ncoeff; ++k) {
+ double A = 1, B = 1;
+ {
+ int compl1 = metricOrder;
+ int compl2 = metricOrder;
+ int compl3 = metricOrder;
+ int compltot = 3*metricOrder;
+ for (int l = 0; l < dimSimplex; l++) {
+ A *= nChoosek(compl1, (int) exp(i, l))
+ * nChoosek(compl2, (int) exp(j, l))
+ * nChoosek(compl3, (int) exp(k, l));
+ B *= nChoosek(compltot, (int) exp(i, l) + (int) exp(j, l) + (int) exp(k, l));
+ compl1 -= (int) exp(i, l);
+ compl2 -= (int) exp(j, l);
+ compl3 -= (int) exp(k, l);
+ compltot -= (int) exp(i, l) + (int) exp(j, l) + (int) exp(k, l);
+ }
+ }
+ for (int l = dimSimplex; l < dim; l++) {
+ A *= nChoosek(metricOrder, (int) exp(i, l))
+ * nChoosek(metricOrder, (int) exp(j, l))
+ * nChoosek(metricOrder, (int) exp(k, l));
+ B *= nChoosek(3*metricOrder, (int) exp(i, l) + (int) exp(j, l) + (int) exp(k, l));
+ }
+
+ if (i == j) {
+ if (j != k) A *= 3;
+ }
+ else {
+ if (j == k || i == k) A *= 3;
+ else A *= 6;
+ }
+
+ ++countP3;
+ int hash = 0;
+ for (int l = 0; l < dim; l++) {
+ hash += (int) (exp(i, l)+exp(j, l)+exp(k, l)) * pow_int(3*metricOrder+1, l);
+ }
+ _ineqP3[hash].push_back(IneqData(A/B, i, j, k));
+ //Msg::Info(" %d", _ineqP3[hash].back().k);
+ }
+ }
+ }
+ _inequations.resize(v.size(), 2);
+ exp = _jacobian->bezier->_exponents;
+ for (unsigned int k = 0; k < v.size(); ++k) {
+ //Msg::Info("%d %d", v[k].first, v[k].second);
+ _inequations(k, 0) = v[k].first;
+ _inequations(k, 1) = v[k].second;
+ }
+
+ int njac = _jacobian->getNumJacNodes();
+ for (int i = 0; i < njac; i++) {
+ for (int j = i; j < njac; j++) {
+ int order = metricOrder/2*3;
+ double A = 1, B = 1;
+ {
+ int compl1 = order;
+ int compl2 = order;
+ int compltot = 2*order;
+ for (int k = 0; k < dimSimplex; k++) {
+ A *= nChoosek(compl1, (int) exp(i, k))
+ * nChoosek(compl2, (int) exp(j, k));
+ B *= nChoosek(compltot, (int) exp(i, k) + (int) exp(j, k));
+ compl1 -= (int) exp(i, k);
+ compl2 -= (int) exp(j, k);
+ compltot -= (int) exp(i, k) + (int) exp(j, k);
+ }
+ }
+ for (int k = dimSimplex; k < dim; k++) {
+ A *= nChoosek(order, (int) exp(i, k))
+ * nChoosek(order, (int) exp(j, k));
+ B *= nChoosek(2*order, (int) exp(i, k) + (int) exp(j, k));
+ }
+
+ if (i != j) A *= 2;
+
+ ++countJ2;
+ int hash = 0;
+ for (int k = 0; k < dim; k++) {
+ hash += (int) (exp(i, k)+exp(j, k)) * pow_int(2*order+1, k);
+ }
+ _ineqJ2[hash].push_back(IneqData(A/B, i, j));
+ }
+ }
+
+ int ret = _lightenInequalities(tol, tol2, countJ2, countP3);
+ if (ret) {
+ Msg::Warning("Heavy Inequalities (error %d)", ret);
+ }
+
+ Msg::Info("A : %d / %d", v.size(), ncoeff*(ncoeff+1)/2);
+ Msg::Info("J2 : %d / %d", countJ2, njac*(njac+1)/2);
+ Msg::Info("P3 : %d / %d", countP3, ncoeff*(ncoeff+1)*(ncoeff+2)/6);
+}
+
+int MetricCoefficient::_lightenInequalities(double tol, double tol2, int &countj, int &countp)
+{
+#ifndef BEFORE
+ 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] = _ineqJ2.begin();
+ itend[0] = _ineqJ2.end();
+ itend[1] = _ineqP3.end();
+ //itend[2] = _ineqJ2.end();
+ for (int k = 0; k < 2; ++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];
+ return 0;
+
+#else
+ std::map<int, std::vector<IneqData> >::iterator itJ, itP;
+ itJ = _ineqJ2.begin();
+ itP = _ineqP3.begin();
+
+ while (itJ != _ineqJ2.end() && itP != _ineqP3.end()) {
+ double verif = 0;
+ for (unsigned int i = 0; i < itJ->second.size(); ++i) {
+ verif += itJ->second[i].val;
+ }
+ if (verif < 1. - 1e-7 || verif > 1. + 1e-7) Msg::Warning("verif J %g", verif);
+ verif = 0;
+ for (unsigned int i = 0; i < itP->second.size(); ++i) {
+ verif += itP->second[i].val;
+ }
+ if (verif < 1.-1e-7 || verif > 1.+1e-7) Msg::Warning("verif P %g", verif);
+ // TODO : remove
+
+ std::sort(itJ->second.begin(), itJ->second.end(), gterIneq());
+ std::sort(itP->second.begin(), itP->second.end(), gterIneq());
+
+ bool print = false;
+ if (true || itJ->second.back().val < tol && itP->second.back().val < tol) {
+ print = true;
+ Msg::Info(" ");
+ //Msg::Info("sizes %d %d", itJ->second.size(), itP->second.size());
+ Msg::Info("JJ-before-JJJJJJJJJJJJJJJJJJJ");
+ for (unsigned int i = 0; i < itJ->second.size(); ++i) {
+ Msg::Info("%g", itJ->second[i].val / tol);
+ }
+ Msg::Info("PP-before-PPPPPPPPPPPPPPPPPPP");
+ for (unsigned int i = 0; i < itP->second.size(); ++i) {
+ Msg::Info("%g", itP->second[i].val / tol);
+ }
+ }
+
+#if defined(HAVE_SCIP)
+ SCIP* scip;
+ SCIP_CALL(SCIPcreate(&scip));
+ SCIP_CALL(SCIPincludeDefaultPlugins(scip));
+ //SCIP_CALL(SCIPsetMessagehdlr(NULL, NULL));
+ SCIP_CALL(SCIPcreateProb(scip, "inequalities", NULL, NULL,
+ NULL, NULL, NULL, NULL, NULL));
+ SCIP_CALL(SCIPsetObjsense(scip, SCIP_OBJSENSE_MAXIMIZE));
+
+ int nj = itJ->second.size();
+ int np = itP->second.size();
+
+ SCIP_VAR** varj = new SCIP_VAR*[nj];
+ SCIP_VAR** varp = new SCIP_VAR*[np];
+ double* cj = new double[nj];
+ double* cp = new double[np];
+
+ for (int i = 0; i < nj; ++i) {
+ std::ostringstream namebuf;
+ namebuf.str("");
+ namebuf << "j" << i;
+
+ SCIP_CALL(SCIPcreateVar(scip, &varj[i], namebuf.str().c_str(), .0, 1.,
+ 1., SCIP_VARTYPE_BINARY, TRUE, FALSE,
+ NULL, NULL, NULL, NULL, NULL));
+ SCIP_CALL(SCIPaddVar(scip, varj[i]));
+ cj[i] = itJ->second[i].val;
+ }
+ for (int i = 0; i < np; ++i) {
+ std::ostringstream namebuf;
+ namebuf.str("");
+ namebuf << "p" << i;
+
+ SCIP_CALL(SCIPcreateVar(scip, &varp[i], namebuf.str().c_str(), .0, 1.,
+ 1., SCIP_VARTYPE_BINARY, TRUE, FALSE,
+ NULL, NULL, NULL, NULL, NULL));
+ SCIP_CALL(SCIPaddVar(scip, varp[i]));
+ cp[i] = itP->second[i].val;
+ }
+
+ SCIP_CONS* cons[3];
+ for (int i = 0; i < nj; ++i) {
+ Msg::Info("+cj%d (%g)", i, cj[i]);
+ }
+
+ Msg::Info("<= tol %g", tol);
+ SCIP_CALL(SCIPcreateConsLinear(scip, &cons[0], "bound J",
+ nj, varj, cj, .0, tol, TRUE,
+ TRUE, TRUE, TRUE, TRUE, FALSE,
+ FALSE, FALSE, FALSE, FALSE));
+ SCIPconsPrint(cons[0], scip->set, NULL, NULL);
+ SCIP_CALL(SCIPcreateConsLinear(scip, &cons[1], "bound P",
+ np, varp, cp, .0, tol, TRUE,
+ TRUE, TRUE, TRUE, TRUE, FALSE,
+ FALSE, FALSE, FALSE, FALSE));
+ SCIP_CALL(SCIPcreateConsLinear(scip, &cons[2], "difference",
+ nj, varj, cj, -tol2, tol2, TRUE,
+ TRUE, TRUE, TRUE, TRUE, FALSE,
+ FALSE, FALSE, FALSE, FALSE));
+ for (int i = 0; i < np; ++i) {
+ //SCIP_CALL(SCIPaddCoefLinear(scip, cons[2], varp[i], -cp[i]));
+ }
+
+ SCIP_CALL(SCIPsolve(scip));
+ SCIPconsPrint(cons[0], scip->set, NULL, NULL);
+ SCIP_SOL* sol = SCIPgetBestSol(scip);
+ if (!sol) {
+ Msg::Error("Scip did not find a solution");
+ }
+ else {
+ int ind = 0;
+ for (int i = 0; i < nj; ++i) {
+ Msg::Info("j%d = %g", i, SCIPgetSolVal(scip, sol, varj[i]));
+ if (SCIPgetSolVal(scip, sol, varj[i]) > .5) {
+ itJ->second[ind] = itJ->second.back();
+ itJ->second.pop_back();
+ }
+ else ++ind;
+ }
+ ind = 0;
+ for (int i = 0; i < np; ++i) {
+ if (SCIPgetSolVal(scip, sol, varp[i]) > .5) {
+ itP->second[ind] = itP->second.back();
+ itP->second.pop_back();
+ }
+ else ++ind;
+ }
+ }
+
+ for (int i = 0; i < nj; ++i)
+ SCIP_CALL(SCIPreleaseVar(scip, &varj[i]));
+ for (int i = 0; i < np; ++i)
+ SCIP_CALL(SCIPreleaseVar(scip, &varp[i]));
+ //SCIP_CALL(SCIPreleaseCons(scip, &cons[0]));
+ //SCIP_CALL(SCIPreleaseCons(scip, &cons[1]));
+ //SCIP_CALL(SCIPreleaseCons(scip, &cons[2]));
+ SCIP_CALL(SCIPfree(&scip));
+ delete[] varj;
+ delete[] varp;
+ delete[] cj;
+ delete[] cp;
+
+#else
+
+ double rmvedJ = 0, rmvedP = 0;
+ bool changed = true;
+
+ while (changed) {
+ changed = false;
+ double potentialJ = rmvedJ, potentialP = rmvedP;
+ int j = itJ->second.size();
+ int p = itP->second.size();
+ bool found = false;
+
+ while (j > 1 && !found && potentialJ <= tol) {
+ double tmp = itP->second[p-1].val;
+ //if (potentialP - potentialJ > -tol2) found = true;
+ while (potentialP + tmp <= tol && potentialP + tmp - potentialJ <= tol2) {
+ potentialP += tmp;
+ --p;
+ if (p < 1) Msg::Fatal("me");
+ tmp = itP->second[p-1].val;
+ if (potentialP - potentialJ >= -tol2) found = true;
+ }
+ if (!found) {
+ potentialJ += itJ->second[j-1].val;
+ --j;
+ if (j < 1) Msg::Fatal("ma");
+ potentialP = rmvedP;
+ p = itP->second.size();
+ }
+ }
+ if (found) {
+ rmvedJ = potentialJ, rmvedP = potentialP;
+ //countJ2 -= itJ->second.size() - j;
+ //countP3 -= itP->second.size() - p; // TODO
+ itJ->second.erase(itJ->second.begin() + j, itJ->second.end());
+ //Msg::Info("%d %d", p, itP->second.size());
+ itP->second.erase(itP->second.begin() + p, itP->second.end());
+ changed = true;
+ }
+ }
+#endif
+
+
+ if (print) {
+ //Msg::Info("sizes %d %d", itJ->second.size(), itP->second.size());
+ Msg::Info("JJ-after-JJJJJJJJJJJJJJJJJJJ");
+ for (unsigned int i = 0; i < itJ->second.size(); ++i) {
+ Msg::Info("%g", itJ->second[i].val / tol);
+ }
+ Msg::Info("PP-after-PPPPPPPPPPPPPPPPPPP");
+ for (unsigned int i = 0; i < itP->second.size(); ++i) {
+ Msg::Info("%g", itP->second[i].val / tol);
+ }
+ }
+
+ ++itJ;
+ ++itP;
+ }
+ return 0;
+#endif
}
void MetricCoefficient::getCoefficients(fullMatrix<double> &coeff, bool bezier)
@@ -112,12 +531,10 @@ void MetricCoefficient::getCoefficients(fullMatrix<double> &coeff, bool bezier)
}
coeff.resize(ref->size1(), 2);
for (int i = 0; i < ref->size1(); ++i) {
- double tmp = pow((*ref)(i, 1), 2);
- tmp += pow((*ref)(i, 2), 2);
- tmp += pow((*ref)(i, 3), 2);
- tmp += 2 * pow((*ref)(i, 4), 2);
- tmp += 2 * pow((*ref)(i, 5), 2);
- tmp += 2 * pow((*ref)(i, 6), 2);
+ double tmp = 0;
+ for (int j = 1; j < 7; ++j) {
+ tmp += pow((*ref)(i, j), 2);
+ }
tmp = std::sqrt(tmp);
double factor = std::sqrt(6)/3;
coeff(i, 0) = (*ref)(i, 0) - factor * tmp;
@@ -131,24 +548,6 @@ void MetricCoefficient::getCoefficients(fullMatrix<double> &coeff, bool bezier)
}
}
-static int nChoosek(int n, int k)
-{
- if (n < k || k < 0) {
- Msg::Error("Wrong argument for combination.");
- 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;
-}
-
void MetricCoefficient::interpolate(const double *uvw, double *minmaxQ)
{
if (minmaxQ == NULL) {
@@ -243,14 +642,12 @@ void MetricCoefficient::interpolate(const double *uvw, double *minmaxQ)
double tmp = pow(terms[1], 2);
tmp += pow(terms[2], 2);
tmp += pow(terms[3], 2);
- tmp += 2 * pow(terms[4], 2);
- tmp += 2 * pow(terms[5], 2);
- tmp += 2 * pow(terms[6], 2);
+ tmp += pow(terms[4], 2);
+ tmp += pow(terms[5], 2);
+ tmp += pow(terms[6], 2);
tmp = std::sqrt(tmp);
double factor = std::sqrt(6)/3;
if (tmp < 1e-3*terms[0]) {
- static int aa = 0;
- Msg::Info("%d", ++aa);
minmaxQ[0] = terms[0] - factor * tmp;
minmaxQ[1] = terms[0] + factor * tmp;
}
@@ -292,68 +689,46 @@ void MetricCoefficient::interpolate(const double *uvw, double *minmaxQ)
delete [] terms;
}
-double MetricCoefficient::getBoundRmin(double tol)
+double MetricCoefficient::getBoundRmin(double tol, int which)
{
+ if (_coefficientsBez.size2() == 0) _computeBezCoeff();
+
fullMatrix<double> nodes(_element->getNumShapeFunctions(),3);
fullVector<double> jac(_jacobian->getNumJacNodes());
_element->getNodesCoord(nodes);
_jacobian->getSignedJacobian(nodes, jac);
- if(jac.size() != _coefficientsBez.size1())
- Msg::Fatal("not same size, jac %d, coeff %d", jac.size(), _coefficientsBez.size1());
-
- // compute RminLag
- double RminLag = 1.;
- const double factor1 = std::sqrt(6) / 3;
- const double factor2 = std::sqrt(6) * 3;
-
- for (int i = 0; i < _bezier->getNumLagCoeff(); ++i) {
- double q = _coefficientsBez(i, 0);
- double tmp = pow_int(_coefficientsBez(i, 1), 2);
- tmp += pow_int(_coefficientsBez(i, 2), 2);
- tmp += pow_int(_coefficientsBez(i, 3), 2);
- tmp += 2 * pow_int(_coefficientsBez(i, 4), 2);
- tmp += 2 * pow_int(_coefficientsBez(i, 5), 2);
- tmp += 2 * pow_int(_coefficientsBez(i, 6), 2);
- tmp = tmp*factor1;
- if (tmp < 1e-3 * q) {
- tmp = (q - tmp) / (q + tmp);
- RminLag = std::min(RminLag, tmp);
- }
- else {
- double phi = jac(i)*jac(i);
- phi -= q * q * q;
- phi += .5 * q * tmp * tmp;
- phi /= tmp * tmp * tmp;
- phi *= factor2;
- if (phi > 1.1 || phi < -1.1) Msg::Warning("phi %g", phi);
- if (phi > 1) phi = 1;
- if (phi < -1) phi = -1;
- phi = std::acos(phi)/3;
- tmp = (q + tmp * std::cos(phi + 2*M_PI/3)) / (q + tmp * std::cos(phi));
- RminLag = std::min(RminLag, tmp);
- }
- }
+ double RminLag, RminBez;
+ if (which == 1)
+ _computeRmin1(_coefficientsBez, jac, RminLag, RminBez, 0);
+ else
+ _computeRmin2(_coefficientsBez, jac, RminLag, RminBez, 0, which);
- // compute RminBez
- /*double RminBez = 1.;
- double minq = minq();
- double minp = minp();
- double maxq = maxq();
- double maxp = maxp();
- if (minq < 1e-3 * maxp) {
- double tmp = (minq - maxp) / (minq + maxp);
- RminLag = std::min(RminBez, tmp);
+ if (RminLag-RminBez < tol) {
+ Msg::Info("RETURNING %g", RminBez);
+ return RminBez;
}
else {
-
- }*/
-
-
- /*compute RminBez
- if RminLag-RminBez < tol : return RminBez
- return subdivide*/
- return RminLag;
+ fullMatrix<double> *copycoeff = new fullMatrix<double>(_coefficientsBez);
+ fullVector<double> *copyjac = new fullVector<double>(jac);
+ MetricData *md = new MetricData(copycoeff, copyjac, RminBez, 0);
+ __numSubdivision = 0;
+ __numSub.resize(20);
+ for (unsigned int i = 0; i < __numSub.size(); ++i) __numSub[i] = 0;
+ __maxdepth = 0;
+ double time = Cpu();
+ double tt = _subdivideForRmin(md, RminLag, tol, which);
+ Msg::Info("> computation time %g", Cpu() - time);
+ Msg::Info("> maxDepth %d", __maxdepth);
+ Msg::Info("> numSubdivision %d", __numSubdivision);
+ int last = __numSub.size();
+ while (--last > 0 && __numSub[last] == 0);
+ for (unsigned int i = 0; i < last+1; ++i) {
+ Msg::Info("> depth %d: %d", i, __numSub[i]);
+ }
+ Msg::Info("RETURNING %g after subdivision", tt);
+ return tt;
+ }
}
int MetricCoefficient::metricOrder(int tag)
@@ -363,9 +738,12 @@ int MetricCoefficient::metricOrder(int tag)
switch (parentType) {
case TYPE_PNT : return 0;
+
case TYPE_LIN : return order;
+
case TYPE_TRI :
case TYPE_TET : return 2*order-2;
+
case TYPE_QUA :
case TYPE_PRI :
case TYPE_HEX :
@@ -435,3 +813,521 @@ void MetricCoefficient::_computeBezCoeff()
_coefficientsLag.size2());
_bezier->matrixLag2Bez.mult(_coefficientsLag, _coefficientsBez);
}
+
+void MetricCoefficient::_computeRmin1(
+ fullMatrix<double> &coeff, fullVector<double> &jac,
+ double &RminLag, double &RminBez, int depth) const
+{
+ RminLag = 1.;
+ static const double factor1 = std::sqrt(6) / 3;
+ static const double factor2 = std::sqrt(6) * 3;
+
+ for (int i = 0; i < _bezier->getNumLagCoeff(); ++i) {
+ double q = coeff(i, 0);
+ double p = pow_int(coeff(i, 1), 2);
+ p += pow_int(coeff(i, 2), 2);
+ p += pow_int(coeff(i, 3), 2);
+ p += pow_int(coeff(i, 4), 2);
+ p += pow_int(coeff(i, 5), 2);
+ p += pow_int(coeff(i, 6), 2);
+ p = std::sqrt(p);
+ if (p < 1e-3 * q) {
+ p *= factor1;
+ RminLag = std::min(RminLag, std::sqrt((q - p) / (q + p)));
+ }
+ else {
+ double phi = q/p;
+ phi = .5 * phi - pow_int(phi,3);
+ phi += jac(i)/p/p*jac(i)/p;
+ phi *= factor2;
+ if (phi > 1.1 || phi < -1.1) {
+ if (!depth) {
+ Msg::Warning("+ phi %g (jac %g, q %g, p %g)", phi, jac(i), q, p);
+ Msg::Info("%g + %g - %g = %g", jac(i)*jac(i)/p/p/p, .5 * q/p, q*q*q/p/p/p, (jac(i)*jac(i)+.5 * q*p*p-q*q*q)/p/p/p);
+ }
+ else if (depth == 1)
+ Msg::Warning("- phi %g @ %d(%d) (jac %g, q %g, p %g)", phi, depth, i, jac(i), q, p);
+ }
+
+ if (phi > 1) phi = 1;
+ if (phi < -1) phi = -1;
+ phi = std::acos(phi)/3;
+ p *= factor1;
+ double tmp = (q + p * std::cos(phi + 2*M_PI/3)) / (q + p * std::cos(phi));
+ if (tmp < 0) Msg::Fatal("1 s normal ? %g", tmp);
+ RminLag = std::min(RminLag, std::sqrt(tmp));
+ }
+ }
+
+ double minq = _minq(coeff);
+ double maxp = _maxp(coeff);
+ if (maxp < 1e-3 * minq) {
+ maxp *= factor1;
+ double tmp = (minq - maxp) / (minq + maxp);
+ RminBez = std::sqrt(tmp);
+ }
+ else {
+ double minp = _minp(coeff);
+ double maxq = _maxq(coeff);
+
+ double minJ2;
+ double maxJ2;
+ _minMaxJacobianSqr(jac, minJ2, maxJ2);
+
+ double minphi = minJ2/pow_int(maxp, 3);
+ minphi -= pow_int(maxq/minp, 3);
+ minphi += .5 * minq/maxp;
+ minphi *= factor2;
+ if (minphi < -1) minphi = -1;
+ if (minphi > 1) minphi = 1;
+ minphi = std::acos(minphi) / 3 + 2*M_PI/3;
+
+ double maxphi = maxJ2/pow_int(minp, 3);
+ maxphi -= pow_int(minq/maxp, 3);
+ maxphi += .5 * maxq/minp;
+ maxphi *= factor2;
+ if (maxphi < -1.) maxphi = -1.;
+ if (maxphi > 1.) maxphi = 1.;
+ maxphi = std::acos(maxphi) / 3;
+
+ double tmp = minq + factor1 * maxp * std::cos(minphi);
+ tmp /= maxq + factor1 * maxp * std::cos(maxphi);
+ if (tmp < 0) RminBez = 0;
+ else RminBez = std::sqrt(tmp);
+ }
+}
+
+void MetricCoefficient::_computeRmin2(
+ fullMatrix<double> &coeff, fullVector<double> &jac,
+ double &RminLag, double &RminBez, int depth, int which) const
+{
+ RminLag = 1.;
+ static const double factor1 = std::sqrt(6) / 3;
+ static const double factor2 = std::sqrt(6) * 3;
+
+ for (int i = 0; i < _bezier->getNumLagCoeff(); ++i) {
+ double q = coeff(i, 0);
+ double p = pow_int(coeff(i, 1), 2);
+ p += pow_int(coeff(i, 2), 2);
+ p += pow_int(coeff(i, 3), 2);
+ p += pow_int(coeff(i, 4), 2);
+ p += pow_int(coeff(i, 5), 2);
+ p += pow_int(coeff(i, 6), 2);
+ p = std::sqrt(p);
+ if (p < 1e-3 * q) {
+ p *= factor1;
+ RminLag = std::min(RminLag, std::sqrt((q - p) / (q + p)));
+ }
+ else {
+ double phi = q/p;
+ phi = .5 * phi - pow_int(phi,3);
+ phi += jac(i)/p/p*jac(i)/p;
+ phi *= factor2;
+ if (phi > 1.1 || phi < -1.1) {
+ if (!depth) {
+ Msg::Warning("+ phi %g (jac %g, q %g, p %g)", phi, jac(i), q, p);
+ Msg::Info("%g + %g - %g = %g", jac(i)*jac(i)/p/p/p, .5 * q/p, q*q*q/p/p/p, (jac(i)*jac(i)+.5 * q*p*p-q*q*q)/p/p/p);
+ }
+ else if (depth == 1)
+ Msg::Warning("- phi %g @ %d(%d) (jac %g, q %g, p %g)", phi, depth, i, jac(i), q, p);
+ }
+
+ if (phi > 1) phi = 1;
+ if (phi < -1) phi = -1;
+ phi = std::acos(phi)/3;
+ p *= factor1;
+ double tmp = (q + p * std::cos(phi + 2*M_PI/3)) / (q + p * std::cos(phi));
+ if (tmp < 0) Msg::Fatal("2 s normal ? %g", tmp);
+ RminLag = std::min(RminLag, std::sqrt(tmp));
+ }
+ }
+
+ double minq = _minq(coeff);
+ double maxp = _maxp(coeff);
+ //Msg::Info("1: %g", maxp/minq);
+ if (maxp < 1e-3 * minq) {
+ maxp *= factor1;
+ double tmp = (minq - maxp) / (minq + maxp); //TODO il y a mieux !
+ RminBez = std::sqrt(tmp);
+ /*if (RminBez > .6 || isnan(RminBez)) Msg::Info("minq %g maxp %g => %g", minq, maxp, RminBez);*/
+ }
+ else {
+ double minJ;
+ if (which == 2 || which == 4) {
+ double maxJ;
+ _minMaxJacobianSqr(jac, minJ, maxJ);
+ minJ /= maxp*maxp*maxp;
+ }
+ else if (which == 3 || which == 5) {
+ _minJ2P3(coeff, jac, minJ);
+ }
+ else {
+ Msg::Fatal("don't know what to do %d", which);
+ return;
+ }
+
+ double a1, a0;
+ {
+ double p = -1./2;
+ double q = -minJ-1/factor2;
+ a1 = cubicCardanoRoot(p, q); //plus grand => -1
+ q = -minJ+1/factor2;
+ a0 = cubicCardanoRoot(p, q); //plus petit => 1
+ }
+
+ double mina;
+ double maxa;
+ if (which == 2 || which == 5) {
+ mina = minq/maxp;
+ double minp = _minp(coeff);
+ if (minp > 0.)
+ maxa = _maxq(coeff)/minp;
+ else
+ maxa = a1;
+ }
+ else if (which == 3 || which == 4) {
+ _minMaxA2(coeff, mina, maxa);
+ }
+ else {
+ Msg::Fatal("don't know what to do %d", which);
+ return;
+ }
+
+ double phimin = std::acos(-factor1/2/mina) - M_PI/3;
+ if (mina > a1) {
+ RminBez = (mina+factor1*std::cos(phimin+2*M_PI/3))/(mina+factor1*std::cos(phimin)); //*
+ RminBez = std::sqrt(RminBez);
+ return;
+ }
+
+ double phim = std::acos(-factor1/2/a1) - M_PI/3;
+ double am;
+ {
+ double p = -1./2;
+ double q = -minJ+std::cos(3*phim)/factor2;
+ am = cubicCardanoRoot(p, q); //milieu
+ }
+
+ if (maxa < am) {
+ double phi = std::acos(factor2*(minJ-maxa*maxa*maxa+maxa/2))/3;
+ RminBez = (maxa+factor1*std::cos(phi+2*M_PI/3))/(maxa+factor1*std::cos(phi));
+ RminBez = std::sqrt(RminBez);
+ return;
+ }
+
+ if (mina > am) {
+ RminBez = (mina+factor1*std::cos(phimin+2*M_PI/3))/(mina+factor1*std::cos(phimin)); //same as *
+ RminBez = std::sqrt(RminBez);
+ return;
+ }
+ else {
+ RminBez = (a1+factor1*std::cos(phim+2*M_PI/3))/(a1+factor1*std::cos(phim));
+ RminBez = std::sqrt(RminBez);
+ return;
+ }
+ }
+}
+
+double MetricCoefficient::_subdivideForRmin(
+ MetricData *md, double RminLag, double tol, int which) const
+{
+ std::priority_queue<MetricData*, std::vector<MetricData*>, lessMinB> subdomains;
+ std::vector<MetricData*> vect;
+ const int numCoeff = md->_metcoeffs->size2();
+ const int numMetPnts = md->_metcoeffs->size1();
+ const int numJacPnts = md->_jaccoeffs->size();
+ const int numSub = _jacobian->getNumDivisions();
+ subdomains.push(md);
+
+ std::vector<fullVector<double>*> trash;
+
+ while (RminLag - subdomains.top()->_RminBez > tol && subdomains.size() < pow_int(8,8)) {
+ fullMatrix<double> *subcoeffs, *coeff;
+ fullVector<double> *subjac, *jac;
+
+ MetricData *current = subdomains.top();
+ subcoeffs = new fullMatrix<double>(numSub*numMetPnts, numCoeff);
+ subjac = new fullVector<double>(numSub*numJacPnts);
+ _bezier->subDivisor.mult(*current->_metcoeffs, *subcoeffs);
+ _jacobian->subdivideBezierCoeff(*current->_jaccoeffs, *subjac);
+ int depth = current->_depth;
+ //Msg::Info("d %d RminBez %g / %g", depth, current->_RminBez, RminLag);
+ delete current;
+ subdomains.pop();
+
+ ++((MetricCoefficient*)this)->__numSubdivision;
+ ++((MetricCoefficient*)this)->__numSub[depth];
+ ((MetricCoefficient*)this)->__maxdepth = std::max(__maxdepth, depth+1);
+
+ if (depth == 0) {
+ Msg::Info("Begin 1");
+ }
+ for (int i = 0; i < numSub; ++i) {
+ coeff = new fullMatrix<double>(numMetPnts, numCoeff);
+ coeff->copy(*subcoeffs, i * numMetPnts, numMetPnts, 0, numCoeff, 0, 0);
+ jac = new fullVector<double>;
+ jac->setAsProxy(*subjac, i * numJacPnts, numJacPnts);
+ double minLag, minBez;
+ if (which == 1)
+ _computeRmin1(*coeff, *jac, minLag, minBez, depth+1);
+ else
+ _computeRmin2(*coeff, *jac, minLag, minBez, depth+1, which);
+ //Msg::Info("new RminBez %g", minBez);
+ RminLag = std::min(RminLag, minLag);
+ MetricData *metData = new MetricData(coeff, jac, minBez, depth+1);
+ subdomains.push(metData);
+ vect.push_back(metData);
+ }
+ if (depth == 0) {
+ Msg::Info("End 1");
+ }
+ trash.push_back(subjac);
+ delete subcoeffs;
+
+ /*for (unsigned int i = 0; i < vect.size(); ++i) {
+ Msg::Info("v %g", vect[i]->_RminBez);
+ }
+ Msg::Info("top %g (RminLag %g)", subdomains.top()->_RminBez, RminLag);
+ return 0;*/
+ //Msg::Info("RminLag %g - RminBez %g @ %d", RminLag, subdomains.top()->_RminBez, subdomains.top()->_depth);
+ }
+
+ md = subdomains.top();
+ double ans = md->_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;
+}
+
+double MetricCoefficient::_minp(const fullMatrix<double> &coeff) const
+{
+ fullMatrix<double> minmaxCoeff(2, 6);
+ for (int j = 0; j < 6; ++j) {
+ minmaxCoeff(0, j) = coeff(0, j+1);
+ minmaxCoeff(1, j) = coeff(0, j+1);
+ }
+
+ for (int i = 1; i < coeff.size1(); ++i) {
+ for (int j = 0; j < 6; ++j) {
+ if (minmaxCoeff(0, j) > coeff(i, j+1)) {
+ minmaxCoeff(0, j) = coeff(i, j+1);
+ }
+ if (minmaxCoeff(1, j) < coeff(i, j+1)) {
+ minmaxCoeff(1, j) = coeff(i, j+1);
+ }
+ }
+ }
+
+ double ans = 0;
+ for (int j = 0; j < 6; ++j) {
+ if (minmaxCoeff(0, j) * minmaxCoeff(1, j) > 0) {
+ ans += minmaxCoeff(0, j) > 0 ?
+ pow_int(minmaxCoeff(0, j), 2) :
+ pow_int(minmaxCoeff(1, j), 2);
+ }
+ }
+ return std::sqrt(ans);
+}
+
+double MetricCoefficient::_minq(const fullMatrix<double> &coeff) const
+{
+ double ans = coeff(0, 0);
+ for (int i = 1; i < coeff.size1(); ++i) {
+ if (ans > coeff(i, 0)) ans = coeff(i, 0);
+ }
+ return ans;
+}
+
+double MetricCoefficient::_maxp(const fullMatrix<double> &coeff) const
+{
+ double ans = 0;
+ for (int i = 0; i < coeff.size1(); ++i) {
+ double tmp = pow_int(coeff(i, 1), 2);
+ tmp += pow_int(coeff(i, 2), 2);
+ tmp += pow_int(coeff(i, 3), 2);
+ tmp += pow_int(coeff(i, 4), 2);
+ tmp += pow_int(coeff(i, 5), 2);
+ tmp += pow_int(coeff(i, 6), 2);
+ if (ans < tmp) ans = tmp;
+ }
+ return std::sqrt(ans);
+}
+
+double MetricCoefficient::_maxq(const fullMatrix<double> &coeff) const
+{
+ double ans = coeff(0, 0);
+ for (int i = 1; i < coeff.size1(); ++i) {
+ if (ans < coeff(i, 0)) ans = coeff(i, 0);
+ }
+ return ans;
+}
+
+void MetricCoefficient::_minMaxA(
+ const fullMatrix<double> &coeff, double &min, double &max) const
+{
+ int i = _inequations(0, 0);
+ int j = _inequations(0, 1);
+
+ // TODO : si coeff(0, m) < 0
+
+ min = max = 0;
+ for (int l = 1; l < 7; ++l) {
+ min += coeff(i, l) * coeff(j, l);
+ }
+ min /= coeff(i, 0) * coeff(j, 0);
+ // TODO : si coeff(0, i|j) == 0
+ max = min;
+
+ for (int k = 1; k < _inequations.size1(); ++k) {
+ double tmp = 0;
+ i = _inequations(k, 0);
+ j = _inequations(k, 1);
+ for (int l = 1; l < 7; ++l) {
+ tmp += coeff(i, l) * coeff(j, l);
+ }
+ tmp /= coeff(i, 0) * coeff(j, 0);
+ min = std::min(tmp, min);
+ max = std::max(tmp, max);
+ }
+
+ if (min < 0)
+ Msg::Fatal("minA ne devrait pas etre negatif");
+
+ min = std::sqrt(min);
+ max = std::sqrt(max);
+ double tmp = min;
+ min = 1./max;
+ max = 1./tmp;
+}
+
+void MetricCoefficient::_minMaxA2(
+ const fullMatrix<double> &coeff, double &min, double &max) const
+{
+ min = 1e10;
+ max = 0;
+ std::map<std::pair<int, int>, double>::const_iterator it;
+
+ for (unsigned int k = 0; k < _ineq2.size(); ++k) {
+ double tmpNum = 0;
+ double tmpDen = 0;
+ //Msg::Info("%d/%d: %d", k, _ineq2.size(), _ineq2[k].size());
+ //Msg::Info(" ");
+ for (it = _ineq2[k].begin(); it != _ineq2[k].end(); ++it) {
+ const int i = it->first.first;
+ const int j = it->first.second;
+ double tmp2 = 0;
+ for (int l = 1; l < 7; ++l) {
+ tmp2 += coeff(i, l) * coeff(j, l);
+ }
+ tmpNum += it->second * tmp2;
+ tmpDen += it->second * coeff(i, 0) * coeff(j, 0);
+ }
+ double val = tmpNum/tmpDen;
+ //Msg::Info("%g/%g = %g", tmpNum, tmpDen, val);
+ min = std::min(val, min);
+ max = std::max(val, max);
+ }
+
+ if (min < 0)
+ Msg::Fatal("minA ne devrait pas etre negatif");
+
+ min = std::sqrt(min);
+ max = std::sqrt(max);
+ double tmp = min;
+ min = 1./max;
+ max = 1./tmp;
+}
+
+void MetricCoefficient::_minMaxJacobianSqr(
+ const fullVector<double> &jac, double &min, double &max) const
+{
+ static int a = 1;
+ if (++a == 1) {
+ for (int i = 1; i < jac.size(); ++i) {
+ Msg::Info("<%g>", jac(i));
+ }
+ }
+ min = max = jac(0);
+ for (int i = 1; i < jac.size(); ++i) {
+ if (min > jac(i)) min = jac(i);
+ if (max < jac(i)) max = jac(i);
+ }
+
+ if (a == 1) {
+ Msg::Info("%g %g", min, max);
+ }
+
+ if (min*max < 0) {
+ max = max > -min ? max*max : min*min;
+ min = 0;
+ }
+ else {
+ if (max > 0) {
+ max = max*max;
+ min = min*min;
+ }
+ else {
+ double tmp = max;
+ max = min*min;
+ min = tmp*tmp;
+ }
+ }
+}
+
+void MetricCoefficient::_minJ2P3(const fullMatrix<double> &coeff,
+ const fullVector<double> &jac, double &min) const
+{
+ fullVector<double> r(coeff.size1());
+ for (int i = 0; i < coeff.size1(); ++i) {
+ r(i) = 0;
+ for (int l = 1; l < 7; ++l) {
+ r(i) += coeff(i, l) * coeff(i, l);
+ }
+ r(i) = std::sqrt(r(i));
+ }
+
+ min = 1e10;
+ std::map<int, std::vector<IneqData> >::const_iterator itJ, itP;
+ itJ = _ineqJ2.begin();
+ itP = _ineqP3.begin();
+
+ //Msg::Warning("sizes %d %d", _ineqJ2.size(), _ineqP3.size());
+ int count = 0;
+ while (itJ != _ineqJ2.end() && itP != _ineqP3.end()) {
+ if (count >= _ineqJ2.size()) Msg::Fatal("aaargh");
+ if (itJ->first != itP->first) Msg::Fatal("not same hash %d %d", itJ->first, itP->first);
+
+ double num = 0;
+ //Msg::Info("sizej %d", itJ->second.size());
+ for (int l = 0; l < itJ->second.size(); ++l) {
+ const int i = itJ->second[l].i;
+ const int j = itJ->second[l].j;
+ num += itJ->second[l].val * jac(i) * jac(j);
+ }
+
+ double den = 0;
+ //Msg::Info("sizep %d", itP->second.size());
+ for (int l = 0; l < itP->second.size(); ++l) {
+ const int i = itP->second[l].i;
+ const int j = itP->second[l].j;
+ const int k = itP->second[l].k;
+ //Msg::Info("i%d j%d k%d", i, j, 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);
+ ++itJ;
+ ++itP;
+ ++count;
+ }
+}
diff --git a/Numeric/MetricBasis.h b/Numeric/MetricBasis.h
index dcf1ab97d919af34e74651c113f720e6627a9a2e..2ea5729fdb38740a7af3bc63e64b3a1985e1824c 100644
--- a/Numeric/MetricBasis.h
+++ b/Numeric/MetricBasis.h
@@ -11,25 +11,82 @@
#include "fullMatrix.h"
class MetricCoefficient {
+private:
+ class MetricData {
+ public:
+ fullMatrix<double> *_metcoeffs;
+ fullVector<double> *_jaccoeffs;
+ double _RminBez;
+ int _depth;
+
+ public:
+ MetricData(fullMatrix<double> *m, fullVector<double> *j, double r, int d) :
+ _metcoeffs(m), _jaccoeffs(j), _RminBez(r), _depth(d) {}
+ ~MetricData() {
+ delete _metcoeffs;
+ delete _jaccoeffs;
+ }
+ };
+ struct lessMinB {
+ bool operator()(const MetricData *md1, const MetricData *md2) const {
+ return md1->_RminBez > md2->_RminBez;
+ }
+ };
+
+ class IneqData {
+ public:
+ int i, j, k;
+ double val;
+
+ public:
+ IneqData(double val, int i, int j, int k = -1) : i(i), j(j), k(k), val(val) {}
+ };
+ class gterIneq {
+ public:
+ bool operator()(const IneqData &id1, const IneqData &id2) const {
+ return id1.val > id2.val;
+ }
+ };
+
private:
MElement *_element;
const JacobianBasis *_jacobian;
const GradientBasis *_gradients;
const bezierBasis *_bezier;
fullMatrix<double> _coefficientsLag, _coefficientsBez;
+ fullMatrix<double> _inequations;
+ std::vector<std::map<std::pair<int, int>, double> > _ineq2; //TODO : pas top
+ std::map<int, std::vector<IneqData> > _ineqJ2, _ineqP3;
+ int __maxdepth, __numSubdivision;
+ std::vector<int> __numSub;
public:
MetricCoefficient(MElement*);
void getCoefficients(fullMatrix<double>&, bool bezier = true);
void interpolate(const double *uvw, double *minmaxQ);
- double getBoundRmin(double tol);
+ double getBoundRmin(double tol, int which);
static int metricOrder(int tag);
private:
void _computeBezCoeff();
+ void _fillInequalities(int order);
+ int _lightenInequalities(double, double, int&, int&); //TODO change
void _interpolateBezierPyramid(const double *uvw, double *minmaxQ);
+ void _computeRmin1(fullMatrix<double>&, fullVector<double>&,
+ double &RminLag, double &RminBez, int) const;
+ void _computeRmin2(fullMatrix<double>&, fullVector<double>&,
+ double &RminLag, double &RminBez, int depth, int which) const;
+ double _subdivideForRmin(MetricData*, double RminLag, double tol, int which) const;
+ double _minp(const fullMatrix<double>&) const;
+ double _minq(const fullMatrix<double>&) const;
+ double _maxp(const fullMatrix<double>&) const;
+ double _maxq(const fullMatrix<double>&) const;
+ void _minMaxA(const fullMatrix<double>&, double &min, double &max) const;
+ void _minMaxA2(const fullMatrix<double>&, double &min, double &max) const;
+ void _minMaxJacobianSqr(const fullVector<double>&, double &min, double &max) const;
+ void _minJ2P3(const fullMatrix<double>&, const fullVector<double>&, double &min) const;
};
#endif
diff --git a/Numeric/bezierBasis.cpp b/Numeric/bezierBasis.cpp
index 3bc3f9c50c6aa2196deb00ac6fda6f9f475299ee..9903be53515a43ea52a1cc4a661537ce164e30e8 100644
--- a/Numeric/bezierBasis.cpp
+++ b/Numeric/bezierBasis.cpp
@@ -551,59 +551,58 @@ void bezierBasis::_construct(int parentType, int p)
return;
}
- int dimSimplex;
switch (parentType) {
case TYPE_PNT :
dim = 0;
numLagCoeff = 1;
- dimSimplex = 0;
+ _dimSimplex = 0;
_exponents = gmshGenerateMonomialsLine(0);
subPoints.push_back(gmshGeneratePointsLine(0));
break;
case TYPE_LIN : {
dim = 1;
- numLagCoeff = 2;
- dimSimplex = 0;
+ numLagCoeff = order == 0 ? 1 : 2;
+ _dimSimplex = 0;
_exponents = gmshGenerateMonomialsLine(order);
subPoints = generateSubPointsLine(order);
break;
}
case TYPE_TRI : {
dim = 2;
- numLagCoeff = 3;
- dimSimplex = 2;
+ numLagCoeff = order == 0 ? 1 : 3;
+ _dimSimplex = 2;
_exponents = gmshGenerateMonomialsTriangle(order);
subPoints = generateSubPointsTriangle(order);
break;
}
case TYPE_QUA : {
dim = 2;
- numLagCoeff = 4;
- dimSimplex = 0;
+ numLagCoeff = order == 0 ? 1 : 4;
+ _dimSimplex = 0;
_exponents = gmshGenerateMonomialsQuadrangle(order);
subPoints = generateSubPointsQuad(order);
break;
}
case TYPE_TET : {
dim = 3;
- numLagCoeff = 4;
- dimSimplex = 3;
+ numLagCoeff = order == 0 ? 1 : 4;
+ _dimSimplex = 3;
_exponents = gmshGenerateMonomialsTetrahedron(order);
subPoints = generateSubPointsTetrahedron(order);
break;
}
case TYPE_PRI : {
dim = 3;
- numLagCoeff = 6;
- dimSimplex = 2;
+ numLagCoeff = order == 0 ? 1 : 6;
+ _dimSimplex = 2;
_exponents = gmshGenerateMonomialsPrism(order);
subPoints = generateSubPointsPrism(order);
break;
}
case TYPE_HEX : {
dim = 3;
- numLagCoeff = 8;
- dimSimplex = 0;
+ numLagCoeff = order == 0 ? 1 : 8;
+ _dimSimplex = 0;
_exponents = gmshGenerateMonomialsHexahedron(order);
subPoints = generateSubPointsHex(order);
break;
@@ -614,7 +613,7 @@ void bezierBasis::_construct(int parentType, int p)
dim = 3;
order = 0;
numLagCoeff = 4;
- dimSimplex = 3;
+ _dimSimplex = 3;
_exponents = gmshGenerateMonomialsTetrahedron(order);
subPoints = generateSubPointsTetrahedron(order);
break;
@@ -625,7 +624,7 @@ void bezierBasis::_construct(int parentType, int p)
fullMatrix<double> bezierPoints = _exponents;
bezierPoints.scale(1./order);
- matrixBez2Lag = generateBez2LagMatrix(_exponents, bezierPoints, order, dimSimplex);
+ matrixBez2Lag = generateBez2LagMatrix(_exponents, bezierPoints, order, _dimSimplex);
matrixBez2Lag.invert(matrixLag2Bez);
- subDivisor = generateSubDivisor(_exponents, subPoints, matrixLag2Bez, order, dimSimplex);
+ subDivisor = generateSubDivisor(_exponents, subPoints, matrixLag2Bez, order, _dimSimplex);
}
diff --git a/Numeric/bezierBasis.h b/Numeric/bezierBasis.h
index 3e236f1aae2ca2a0724bc8d35bec42686b533d29..257e54be4a66424454486c269686f51ba02effd6 100644
--- a/Numeric/bezierBasis.h
+++ b/Numeric/bezierBasis.h
@@ -19,6 +19,8 @@ class bezierBasis {
int numLagCoeff;
int dim, type, order;
int numDivisions;
+ public:
+ int _dimSimplex;
fullMatrix<double> _exponents;
public :