From a590fefd20889fdaad300fbdb1aa895be1cf2306 Mon Sep 17 00:00:00 2001
From: Amaury Johnan <amjohnen@gmail.com>
Date: Mon, 26 May 2014 13:54:03 +0000
Subject: [PATCH] cleanup
---
Numeric/MetricBasis.cpp | 642 +---------------------------------------
Numeric/MetricBasis.h | 51 +---
2 files changed, 6 insertions(+), 687 deletions(-)
diff --git a/Numeric/MetricBasis.cpp b/Numeric/MetricBasis.cpp
index eadcbe3d17..fd15fc294c 100644
--- a/Numeric/MetricBasis.cpp
+++ b/Numeric/MetricBasis.cpp
@@ -172,9 +172,9 @@ double MetricBasis::getMinR(const MElement *el, MetricData *&md, int deg) const
_maxKstA3(*coeff, *jac, mina, beta, maxK3);
if (md->_num == 22)
- _computeRmin3(*coeff, *jac, RminLag, RminBez, 0, true);
+ _computeRmin(*coeff, *jac, RminLag, RminBez, 0, true);
else
- _computeRmin3(*coeff, *jac, RminLag, RminBez, 0, false);
+ _computeRmin(*coeff, *jac, RminLag, RminBez, 0, false);
((MetricBasis*)this)->file << minK*6*std::sqrt(6) << " ";
((MetricBasis*)this)->file << maxJ2/minpp/minpp/minpp*6*std::sqrt(6) << " ";
@@ -290,17 +290,7 @@ double MetricBasis::getBoundRmin(const MElement *el, MetricData *&md, fullMatrix
}
Msg::Info("----------------");*/
- if (MetricBasis::_which == 1)
- _computeRmin1(*metCoeff, *jac, RminLag, RminBez, 0);
- else if (MetricBasis::_which == 0) {
- _computeRmin3(*metCoeff, *jac, RminLag, RminBez, 0);
- /*double tmpL, tmpB;
- _computeRmin2(*metCoeff, *jac, tmpL, tmpB, 0, 3);
- if (std::abs(tmpL-RminLag) > 1e-7) Msg::Warning("Lag %g %g", tmpL, RminLag);
- if (std::abs(tmpB-RminBez) > 1e-3) Msg::Warning("Bez %g %g %g", tmpB, RminBez, tmpL);*/
- }
- else
- _computeRmin2(*metCoeff, *jac, RminLag, RminBez, 0, MetricBasis::_which);
+ _computeRmin(*metCoeff, *jac, RminLag, RminBez, 0);
fullVector<double> *jjac = new fullVector<double>(*jac);
fullMatrix<double> *mmet = new fullMatrix<double>(*metCoeff);
@@ -347,70 +337,6 @@ double MetricBasis::getBoundRmin(const MElement *el, MetricData *&md, fullMatrix
}
}
-MetricCoefficient::MetricCoefficient(MElement *el) : _element(el)
-{
- const int tag = el->getTypeForMSH();
- const int type = ElementType::ParentTypeFromTag(tag);
- const int metricOrder = MetricBasis::metricOrder(tag);
- if (type == TYPE_HEX || type == TYPE_PRI) {
- int order = ElementType::OrderFromTag(tag);
- _jacobian = new JacobianBasis(tag, 3*order);
- }
- 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);
-
- int nSampPnts = _gradients->getNumSamplingPoints();
- int nMapping = _gradients->getNumMapNodes();
- fullMatrix<double> nodesXYZ(nMapping, 3);
- el->getNodesCoord(nodesXYZ);
-
- switch (el->getDim()) {
- case 0 :
- return;
-
- case 1 :
- {
- Msg::Fatal("not implemented");
- }
- break;
-
- case 2 :
- {
- Msg::Fatal("not implemented");
- }
- break;
-
- case 3 :
- {
- fullMatrix<double> dxyzdX(nSampPnts,3), dxyzdY(nSampPnts,3), dxyzdZ(nSampPnts,3);
- _gradients->getGradientsFromNodes(nodesXYZ, &dxyzdX, &dxyzdY, &dxyzdZ);
-
- _coefficientsLag.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;
- _coefficientsLag(i, 0) = (dvxdX + dvxdY + dvxdZ) / 3;
- _coefficientsLag(i, 1) = dvxdX - _coefficientsLag(i, 0);
- _coefficientsLag(i, 2) = dvxdY - _coefficientsLag(i, 0);
- _coefficientsLag(i, 3) = dvxdZ - _coefficientsLag(i, 0);
- 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;
- }
-}
-
void MetricBasis::_fillInequalities(int metricOrder)
{
int dimSimplex = _bezier->_dimSimplex;
@@ -552,7 +478,6 @@ void MetricBasis::_fillInequalities(int metricOrder)
Msg::Info("P3 : %d / %d", countP3, ncoeff*(ncoeff+1)*(ncoeff+2)/6);
}
-
void MetricBasis::_lightenInequalities(int &countj, int &countp, int &counta)
{
double tol = .0;
@@ -588,59 +513,6 @@ void MetricBasis::_lightenInequalities(int &countj, int &countp, int &counta)
counta -= cnt[2];
}
-void MetricCoefficient::getCoefficients(fullMatrix<double> &coeff, bool bezier)
-{
- fullMatrix<double> *ref;
- switch (_coefficientsLag.size2()) {
- case 1:
- if (!bezier) coeff = _coefficientsLag;
- else {
- if (_coefficientsBez.size2() != 1) _computeBezCoeff();
- coeff = _coefficientsBez;
- }
- break;
-
- case 3:
- if (!bezier) ref = &_coefficientsLag;
- else {
- if (_coefficientsBez.size2() != 3) _computeBezCoeff();
- ref = &_coefficientsBez;
- }
- 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 = std::sqrt(tmp);
- coeff(i, 0) = (*ref)(i, 0) - tmp;
- coeff(i, 1) = (*ref)(i, 0) + tmp;
- }
- break;
-
- case 7:
- if (!bezier) ref = &_coefficientsLag;
- else {
- if (_coefficientsBez.size2() != 7) _computeBezCoeff();
- ref = &_coefficientsBez;
- }
- coeff.resize(ref->size1(), 2);
- for (int i = 0; i < ref->size1(); ++i) {
- 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;
- coeff(i, 1) = (*ref)(i, 0) + factor * tmp;
- }
- break;
-
- default:
- Msg::Error("Wrong number of functions for metric: %d",
- _coefficientsLag.size2());
- }
-}
-
void MetricBasis::interpolate(const MElement *el, const MetricData *md, const double *uvw, double *minmaxQ, bool write) const
{
if (minmaxQ == NULL) {
@@ -779,190 +651,6 @@ void MetricBasis::interpolate(const MElement *el, const MetricData *md, const do
delete [] terms;
}
-void MetricCoefficient::interpolate(const double *uvw, double *minmaxQ)
-{
- if (minmaxQ == NULL) {
- Msg::Error("Cannot write solution of interpolation");
- return;
- }
-
- int order = _bezier->getOrder();
-
- int dimSimplex = 0;
- fullMatrix<double> exponents;
- double bezuvw[3];
- switch (_element->getType()) {
- case TYPE_PYR:
- bezuvw[0] = .5 * (1 + uvw[0]);
- bezuvw[1] = .5 * (1 + uvw[1]);
- bezuvw[2] = uvw[2];
- _interpolateBezierPyramid(uvw, minmaxQ);
- return;
-
- case TYPE_HEX:
- bezuvw[0] = .5 * (1 + uvw[0]);
- bezuvw[1] = .5 * (1 + uvw[1]);
- bezuvw[2] = .5 * (1 + uvw[2]);
- dimSimplex = 0;
- exponents = gmshGenerateMonomialsHexahedron(order);
- break;
-
- case TYPE_TET:
- bezuvw[0] = uvw[0];
- bezuvw[1] = uvw[1];
- bezuvw[2] = uvw[2];
- dimSimplex = 3;
- exponents = gmshGenerateMonomialsTetrahedron(order);
- break;
-
- case TYPE_PRI:
- bezuvw[0] = uvw[0];
- bezuvw[1] = uvw[1];
- bezuvw[2] = .5 * (1 + uvw[2]);
- dimSimplex = 2;
- exponents = gmshGenerateMonomialsPrism(order);
- break;
- }
-
- int numCoeff = exponents.size1();
- int dim = exponents.size2();
-
- if (_coefficientsBez.size2() == 0) _computeBezCoeff();
-
- double *terms = new double[_coefficientsBez.size2()];
- for (int t = 0; t < _coefficientsBez.size2(); ++t) {
- terms[t] = 0;
- for (int i = 0; i < numCoeff; i++) {
- double dd = 1;
- double pointCompl = 1.;
- int exponentCompl = order;
- for (int k = 0; k < dimSimplex; k++) {
- dd *= nChoosek(exponentCompl, (int) exponents(i, k))
- * pow(bezuvw[k], exponents(i, k));
- pointCompl -= bezuvw[k];
- exponentCompl -= (int) exponents(i, k);
- }
- dd *= pow(pointCompl, exponentCompl);
-
- for (int k = dimSimplex; k < dim; k++)
- dd *= nChoosek(order, (int) exponents(i, k))
- * pow(bezuvw[k], exponents(i, k))
- * pow(1. - bezuvw[k], order - exponents(i, k));
- terms[t] += _coefficientsBez(i, t) * dd;
- }
- }
-
- switch (_coefficientsBez.size2()) {
- case 1:
- minmaxQ[0] = terms[0];
- minmaxQ[1] = terms[0];
- break;
-
- case 3:
- {
- double tmp = pow(terms[1], 2);
- tmp += pow(terms[2], 2);
- tmp = std::sqrt(tmp);
- minmaxQ[0] = terms[0] - tmp;
- minmaxQ[1] = terms[0] + tmp;
- }
- break;
-
- case 7:
- {
- double tmp = pow(terms[1], 2);
- tmp += pow(terms[2], 2);
- tmp += pow(terms[3], 2);
- tmp += pow(terms[4], 2);
- tmp += pow(terms[5], 2);
- tmp += pow(terms[6], 2);
- tmp = std::sqrt(tmp);
- double factor = std::sqrt(6)/3;
- if (tmp < 1e-3*terms[0]) {
- minmaxQ[0] = terms[0] - factor * tmp;
- minmaxQ[1] = terms[0] + factor * tmp;
- }
- else {
- double phi;
- {
- fullMatrix<double> nodes(_element->getNumShapeFunctions(),3);
- _element->getNodesCoord(nodes);
- fullVector<double> jac(_jacobian->getNumJacNodes());
- _jacobian->getSignedJacobian(nodes, jac);
-
- nodes.resize(1, 3);
- nodes(0, 0) = uvw[0];
- nodes(0, 1) = uvw[1];
- nodes(0, 2) = uvw[2];
-
- fullMatrix<double> result;
- _jacobian->interpolate(jac, nodes, result, false);
- phi = result(0, 0)*result(0, 0);
- }
- phi -= terms[0]*terms[0]*terms[0];
- phi += .5*terms[0]*tmp*tmp;
- phi /= tmp*tmp*tmp;
- phi *= 3*std::sqrt(6);
- if (phi > 1) phi = 1;
- if (phi < -1) phi = -1;
- phi = std::acos(phi)/3;
- minmaxQ[0] = terms[0] + factor * tmp * std::cos(phi + 2*M_PI/3);
- minmaxQ[1] = terms[0] + factor * tmp * std::cos(phi);
- }
- }
- break;
-
- default:
- Msg::Error("Wrong number of functions for metric: %d",
- _coefficientsLag.size2());
- }
-
- delete [] terms;
-}
-
-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);
-
- double RminLag, RminBez;
- if (which == 1)
- _basis->_computeRmin1(_coefficientsBez, jac, RminLag, RminBez, 0);
- else
- _basis->_computeRmin2(_coefficientsBez, jac, RminLag, RminBez, 0, which);
-
- if (RminLag-RminBez < tol) {
- Msg::Info("RETURNING %g", RminBez);
- return RminBez;
- }
- else {
- fullMatrix<double> *copycoeff = new fullMatrix<double>(_coefficientsBez);
- fullVector<double> *copyjac = new fullVector<double>(jac);
- MetricBasis::MetricData *md = new MetricBasis::MetricData(copycoeff, copyjac, RminBez, 0, 0);
- //Msg::Info("+3 %d", md);
- __numSubdivision = 0;
- __numSub.resize(20);
- for (unsigned int i = 0; i < __numSub.size(); ++i) __numSub[i] = 0;
- __maxdepth = 0;
- double time = Cpu();
- double tt = _basis->_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 (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 MetricBasis::metricOrder(int tag)
{
const int parentType = ElementType::ParentTypeFromTag(tag);
@@ -986,322 +674,7 @@ int MetricBasis::metricOrder(int tag)
}
}
-void MetricCoefficient::_interpolateBezierPyramid(const double *uvw, double *minmaxQ)
-{
- int order = _jacobian->bezier->getOrder();
-
- fullMatrix<double> exponents = JacobianBasis::generateJacMonomialsPyramid(order);
- int numCoeff = exponents.size1();
-
- if (_coefficientsBez.size2() == 0) _computeBezCoeff();
-
- static int aa = 0;
- if (++aa == 1) {
- fullMatrix<double> val;
- getCoefficients(val, true);
- for (int i = 0; i < val.size2(); ++i) {
- Msg::Info("--------- column %d", i);
- for (int j = 0; j < val.size1(); ++j) {
- Msg::Info("%.2e", val(j, i));
- }
- }
- }
- double terms[7] = {0,0,0,0,0,0,0};
- for (int t = 0; t < _coefficientsBez.size2(); ++t) {
- terms[t] = 0;
- for (int j = 0; j < numCoeff; j++) {
- double dd =
- nChoosek(order, exponents(j, 2))
- * pow(uvw[2], exponents(j, 2))
- * pow(1. - uvw[2], order - exponents(j, 2));
-
- double p = order + 2 - exponents(j, 2);
- double denom = 1. - uvw[2];
- dd *= nChoosek(p, exponents(j, 0))
- * nChoosek(p, exponents(j, 1))
- * pow(uvw[0] / denom, exponents(j, 0))
- * pow(uvw[1] / denom, exponents(j, 1))
- * pow(1. - uvw[0] / denom, p - exponents(j, 0))
- * pow(1. - uvw[1] / denom, p - exponents(j, 1));
- terms[t] += _coefficientsBez(j, t) * dd;
- }
- }
-
- 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 = std::sqrt(tmp);
- double factor = std::sqrt(6)/3;
- minmaxQ[0] = terms[0] - factor * tmp;
- minmaxQ[1] = terms[0] + factor * tmp;
-}
-
-void MetricCoefficient::_computeBezCoeff()
-{
- _coefficientsBez.resize(_coefficientsLag.size1(),
- _coefficientsLag.size2());
- _bezier->matrixLag2Bez.mult(_coefficientsLag, _coefficientsBez);
-}
-
-void MetricBasis::_computeRmin1(
- const fullMatrix<double> &coeff, const 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 MetricBasis::_computeRmin2(
- const fullMatrix<double> &coeff, const 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) {
- if (tmp < -1e-7) Msg::Fatal("2 s normal ? %g", tmp);
- else tmp = 0;
- }
- 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, maxJ;
- if (which == 2 || which == 4) {
- _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; //TODO : a0 pas n残essaire
- {
- 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, maxa2;
- double minp;
- 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);
- _maxAstK(coeff, jac, minJ, a1, maxa2);
- maxa = maxa2;
- }
- else {
- Msg::Fatal("don't know what to do %d", which);
- 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
- }
-
-
- {
- double pCorner[4] = {0, 0, 0, 0};
- for (int k = 0; k < 4; ++k) {
- for (int i = 1; i < 7; ++i) {
- pCorner[k] += pow_int(coeff(k, i), 2);
- }
- pCorner[k] = std::sqrt(pCorner[k]);
- }
- double lagminJ, lagmina, lagmaxa;
- lagminJ = jac(0)*jac(0) / pCorner[0]/pCorner[0]/pCorner[0];
- lagmina = lagmaxa = coeff(0, 0) / pCorner[0];
- static int aaa = 0;
- //if (aaa == 0) Msg::Info(" J%g a%g", lagminJ, lagmina);
- for (int k = 1; k < 4; ++k) {
- lagminJ = std::min(lagminJ, jac(k)*jac(k) / pCorner[k]/pCorner[k]/pCorner[k]);
- lagmina = std::min(lagmina, coeff(k, 0)/pCorner[k]);
- lagmaxa = std::max(lagmaxa, coeff(k, 0)/pCorner[k]);
- //if (aaa == 0) Msg::Info(" J%g a%g", jac(k)*jac(k) / pCorner[k]/pCorner[k]/pCorner[k], coeff(k, 0)/pCorner[k]);
- }
- ++aaa;
- //Msg::Info("minJ %g[%g] (a0,am,a1)(%g, %g,%g) a(%g[%g],%g[%g],(%g)) beta%g",
- // minJ, lagminJ, a0, am, a1, mina, lagmina, maxa, lagmaxa, maxa2, 1. / (1 - 1./6/a1/a1));
- }
-
-
- if (maxa < am) {
- double phi = std::acos(factor2*(minJ-maxa*maxa*maxa+maxa/2))/3;
- RminBez = (maxa+factor1*std::cos(phi+2*M_PI/3))/(maxa+factor1*std::cos(phi));
- RminBez = std::sqrt(RminBez);
- //Msg::Info("2");
- return;
- }
-
- double phimin = std::acos(-factor1/2/mina) - M_PI/3;
- if (mina > am) {
- if (which != 2 && which != 5) {
- minp = _minp(coeff);
- }
- if (which != 2 && which != 4) {
- double tmp;
- _minMaxJacobianSqr(jac, tmp, maxJ);
- maxJ /= minp*minp*minp;
- }
- double myphi = std::acos(factor2*(maxJ-mina*mina*mina+mina/2))/3;
- myphi = std::max(phimin, myphi);
-
- RminBez = (mina+factor1*std::cos(myphi+2*M_PI/3))/(mina+factor1*std::cos(myphi));
- // TODO verif plus haut si faut pas faire la m仁e chose du coup.
- RminBez = std::sqrt(RminBez);
- //Msg::Info("3 minJ %g maxJ %g am %g a1 %g mina %g maxa %g", minJ, maxJ, am, a1, mina, maxa);
- //Msg::Info("3 phi %g %g %g", phimin, std::acos(factor2*(maxJ-mina*mina*mina+mina/2))/3, std::acos(-1)/3);
- return;
- }
- else {
- RminBez = (am+factor1*std::cos(phim+2*M_PI/3))/(am+factor1*std::cos(phim));
- RminBez = std::sqrt(RminBez);
- //Msg::Info("4");
- return;
- }
- }
-}
-
-void MetricBasis::_computeRmin3(
+void MetricBasis::_computeRmin(
const fullMatrix<double> &coeff, const fullVector<double> &jac,
double &RminLag, double &RminBez,
int depth, bool debug) const
@@ -1512,12 +885,7 @@ double MetricBasis::_subdivideForRmin(
jac = new fullVector<double>;
jac->setAsProxy(*subjac, i * numJacPnts, numJacPnts);
double minLag, minBez;
- if (which == 1)
- _computeRmin1(*coeff, *jac, minLag, minBez, depth+1);
- else if (which == 0)
- _computeRmin3(*coeff, *jac, minLag, minBez, depth+1);
- else
- _computeRmin2(*coeff, *jac, minLag, minBez, depth+1, which);
+ _computeRmin(*coeff, *jac, minLag, minBez, depth+1);
//Msg::Info("new RminBez %g", minBez);
RminLag = std::min(RminLag, minLag);
int newNum = num + (i+1) * pow_int(10, depth);
diff --git a/Numeric/MetricBasis.h b/Numeric/MetricBasis.h
index b57effcf88..dbcc82285d 100644
--- a/Numeric/MetricBasis.h
+++ b/Numeric/MetricBasis.h
@@ -75,11 +75,7 @@ private:
void _fillInequalities(int order);
void _lightenInequalities(int&, int&, int&); //TODO change
- void _computeRmin1(const fullMatrix<double>&, const fullVector<double>&,
- double &RminLag, double &RminBez, int) const;
- void _computeRmin2(const fullMatrix<double>&, const fullVector<double>&,
- double &RminLag, double &RminBez, int depth, int which) const;
- void _computeRmin3(const fullMatrix<double>&, const fullVector<double>&,
+ void _computeRmin(const fullMatrix<double>&, const fullVector<double>&,
double &RminLag, double &RminBez, int depth, bool debug = false) const;
double _subdivideForRmin(MetricData*, double RminLag, double tol, int which) const;
@@ -121,49 +117,4 @@ private:
};
};
-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;
- }
- };
-
- private:
- MElement *_element;
- const JacobianBasis *_jacobian;
- const GradientBasis *_gradients;
- const bezierBasis *_bezier;
- fullMatrix<double> _coefficientsLag, _coefficientsBez;
- int __maxdepth, __numSubdivision;
- std::vector<int> __numSub;
- MetricBasis *_basis;
-
- public:
- MetricCoefficient(MElement*);
-
- void getCoefficients(fullMatrix<double>&, bool bezier = true);
- void interpolate(const double *uvw, double *minmaxQ);
- double getBoundRmin(double tol, int which);
-
- private:
- void _computeBezCoeff();
- void _interpolateBezierPyramid(const double *uvw, double *minmaxQ);
-};
-
#endif
--
GitLab