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

cleanup

parent f7f0fa2c
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Numeric/MetricBasis.cpp
+ 5
− 637
View file @ a590fefd
...@@ -172,9 +172,9 @@ double MetricBasis::getMinR(const MElement *el, MetricData *&md, int deg) const ...@@ -172,9 +172,9 @@ double MetricBasis::getMinR(const MElement *el, MetricData *&md, int deg) const
_maxKstA3(*coeff, *jac, mina, beta, maxK3); _maxKstA3(*coeff, *jac, mina, beta, maxK3);
if (md->_num == 22) if (md->_num == 22)
_computeRmin3(*coeff, *jac, RminLag, RminBez, 0, true); _computeRmin(*coeff, *jac, RminLag, RminBez, 0, true);
else 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 << minK*6*std::sqrt(6) << " ";
((MetricBasis*)this)->file << maxJ2/minpp/minpp/minpp*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 ...@@ -290,17 +290,7 @@ double MetricBasis::getBoundRmin(const MElement *el, MetricData *&md, fullMatrix
} }
Msg::Info("----------------");*/ Msg::Info("----------------");*/
if (MetricBasis::_which == 1) _computeRmin(*metCoeff, *jac, RminLag, RminBez, 0);
_computeRmin1(*metCoeff, *jac, RminLag, RminBez, 0);
else if (MetricBasis::_which == 0) {
_computeRmin3(*metCoeff, *jac, RminLag, RminBez, 0);
/*double tmpL, tmpB;
_computeRmin2(*metCoeff, *jac, tmpL, tmpB, 0, 3);
if (std::abs(tmpL-RminLag) > 1e-7) Msg::Warning("Lag %g %g", tmpL, RminLag);
if (std::abs(tmpB-RminBez) > 1e-3) Msg::Warning("Bez %g %g %g", tmpB, RminBez, tmpL);*/
}
else
_computeRmin2(*metCoeff, *jac, RminLag, RminBez, 0, MetricBasis::_which);
fullVector<double> *jjac = new fullVector<double>(*jac); fullVector<double> *jjac = new fullVector<double>(*jac);
fullMatrix<double> *mmet = new fullMatrix<double>(*metCoeff); fullMatrix<double> *mmet = new fullMatrix<double>(*metCoeff);
...@@ -347,70 +337,6 @@ double MetricBasis::getBoundRmin(const MElement *el, MetricData *&md, fullMatrix ...@@ -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) void MetricBasis::_fillInequalities(int metricOrder)
{ {
int dimSimplex = _bezier->_dimSimplex; int dimSimplex = _bezier->_dimSimplex;
...@@ -552,7 +478,6 @@ void MetricBasis::_fillInequalities(int metricOrder) ...@@ -552,7 +478,6 @@ void MetricBasis::_fillInequalities(int metricOrder)
Msg::Info("P3 : %d / %d", countP3, ncoeff*(ncoeff+1)*(ncoeff+2)/6); Msg::Info("P3 : %d / %d", countP3, ncoeff*(ncoeff+1)*(ncoeff+2)/6);
} }
void MetricBasis::_lightenInequalities(int &countj, int &countp, int &counta) void MetricBasis::_lightenInequalities(int &countj, int &countp, int &counta)
{ {
double tol = .0; double tol = .0;
...@@ -588,59 +513,6 @@ void MetricBasis::_lightenInequalities(int &countj, int &countp, int &counta) ...@@ -588,59 +513,6 @@ void MetricBasis::_lightenInequalities(int &countj, int &countp, int &counta)
counta -= cnt[2]; 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 void MetricBasis::interpolate(const MElement *el, const MetricData *md, const double *uvw, double *minmaxQ, bool write) const
{ {
if (minmaxQ == NULL) { if (minmaxQ == NULL) {
...@@ -779,190 +651,6 @@ void MetricBasis::interpolate(const MElement *el, const MetricData *md, const do ...@@ -779,190 +651,6 @@ void MetricBasis::interpolate(const MElement *el, const MetricData *md, const do
delete [] terms; 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) int MetricBasis::metricOrder(int tag)
{ {
const int parentType = ElementType::ParentTypeFromTag(tag); const int parentType = ElementType::ParentTypeFromTag(tag);
...@@ -986,322 +674,7 @@ int MetricBasis::metricOrder(int tag) ...@@ -986,322 +674,7 @@ int MetricBasis::metricOrder(int tag)
} }
} }
void MetricCoefficient::_interpolateBezierPyramid(const double *uvw, double *minmaxQ) void MetricBasis::_computeRmin(
{
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 ncessaire
{
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 mme chose du coup.
RminBez = std::sqrt(RminBez);
//Msg::Info("3 minJ %g maxJ %g am %g a1 %g mina %g maxa %g", minJ, maxJ, am, a1, mina, maxa);
//Msg::Info("3 phi %g %g %g", phimin, std::acos(factor2*(maxJ-mina*mina*mina+mina/2))/3, std::acos(-1)/3);
return;
}
else {
RminBez = (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(
const fullMatrix<double> &coeff, const fullVector<double> &jac, const fullMatrix<double> &coeff, const fullVector<double> &jac,
double &RminLag, double &RminBez, double &RminLag, double &RminBez,
int depth, bool debug) const int depth, bool debug) const
...@@ -1512,12 +885,7 @@ double MetricBasis::_subdivideForRmin( ...@@ -1512,12 +885,7 @@ double MetricBasis::_subdivideForRmin(
jac = new fullVector<double>; jac = new fullVector<double>;
jac->setAsProxy(*subjac, i * numJacPnts, numJacPnts); jac->setAsProxy(*subjac, i * numJacPnts, numJacPnts);
double minLag, minBez; double minLag, minBez;
if (which == 1) _computeRmin(*coeff, *jac, minLag, minBez, depth+1);
_computeRmin1(*coeff, *jac, minLag, minBez, depth+1);
else if (which == 0)
_computeRmin3(*coeff, *jac, minLag, minBez, depth+1);
else
_computeRmin2(*coeff, *jac, minLag, minBez, depth+1, which);
//Msg::Info("new RminBez %g", minBez); //Msg::Info("new RminBez %g", minBez);
RminLag = std::min(RminLag, minLag); RminLag = std::min(RminLag, minLag);
int newNum = num + (i+1) * pow_int(10, depth); int newNum = num + (i+1) * pow_int(10, depth);
......
Numeric/MetricBasis.h
+ 1
− 50
View file @ a590fefd
...@@ -75,11 +75,7 @@ private: ...@@ -75,11 +75,7 @@ private:
void _fillInequalities(int order); void _fillInequalities(int order);
void _lightenInequalities(int&, int&, int&); //TODO change void _lightenInequalities(int&, int&, int&); //TODO change
void _computeRmin1(const fullMatrix<double>&, const fullVector<double>&, void _computeRmin(const fullMatrix<double>&, const fullVector<double>&,
double &RminLag, double &RminBez, int) const;
void _computeRmin2(const fullMatrix<double>&, const fullVector<double>&,
double &RminLag, double &RminBez, int depth, int which) const;
void _computeRmin3(const fullMatrix<double>&, const fullVector<double>&,
double &RminLag, double &RminBez, int depth, bool debug = false) const; double &RminLag, double &RminBez, int depth, bool debug = false) const;
double _subdivideForRmin(MetricData*, double RminLag, double tol, int which) const; double _subdivideForRmin(MetricData*, double RminLag, double tol, int which) const;
...@@ -121,49 +117,4 @@ private: ...@@ -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 #endif
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment