diff --git a/Numeric/JacobianBasis.cpp b/Numeric/JacobianBasis.cpp
index f39583a4e532ee971dcc97880682dfab16e8aa02..13320799a460b31a27b10ec0e42a206d684f0434 100644
--- a/Numeric/JacobianBasis.cpp
+++ b/Numeric/JacobianBasis.cpp
@@ -25,7 +25,6 @@ namespace {
GradientBasis::GradientBasis(int tag, int order)
{
- _bezier = BasisFactory::getBezierBasis(ElementType::ParentTypeFromTag(tag), order);
const int type = ElementType::ParentTypeFromTag(tag);
fullMatrix<double> samplingPoints;
@@ -57,7 +56,7 @@ GradientBasis::GradientBasis(int tag, int order)
//samplingPoints = JacobianBasis::generateJacPointsPyramid(order);
break;
default :
- Msg::Error("Unknown Jacobian function space for element type %d", type);
+ Msg::Error("Unknown Jacobian function space for element tag %d", tag);
return;
}
const int numSampPnts = samplingPoints.size1();
@@ -125,7 +124,7 @@ JacobianBasis::JacobianBasis(int tag)
lagPoints = generateJacPointsPyramid(jacobianOrder);
break;
default :
- Msg::Error("Unknown Jacobian function space for element type %d", parentType);
+ Msg::Error("Unknown Jacobian function space for element tag %d", tag);
return;
}
numJacNodes = lagPoints.size1();
@@ -134,7 +133,7 @@ JacobianBasis::JacobianBasis(int tag)
bezier = BasisFactory::getBezierBasis(parentType, jacobianOrder);
// Store shape function gradients of mapping at Jacobian nodes
- _gradBasis = new GradientBasis(tag, jacobianOrder);
+ _gradBasis = BasisFactory::getGradientBasis(tag, jacobianOrder);
// Compute matrix for lifting from primary Jacobian basis to Jacobian basis
int primJacType = ElementType::getTag(parentType, primJacobianOrder, false);
@@ -626,8 +625,8 @@ void JacobianBasis::getMetricMinAndGradients(const fullMatrix<double> &nodesXYZ,
for (int l = 0; l < numJacNodes; l++) {
double jac[2][2] = {{0., 0.}, {0., 0.}};
for (int i = 0; i < numMapNodes; i++) {
- const double &dPhidX = _gradBasis->getgradShapeMatX(l,i);
- const double &dPhidY = _gradBasis->getgradShapeMatY(l,i);
+ const double &dPhidX = _gradBasis->gradShapeMatX(l,i);
+ const double &dPhidY = _gradBasis->gradShapeMatY(l,i);
const double dpsidx = dPhidX * invJaci[0][0] + dPhidY * invJaci[1][0];
const double dpsidy = dPhidX * invJaci[0][1] + dPhidY * invJaci[1][1];
jac[0][0] += nodesXYZ(i,0) * dpsidx;
diff --git a/Numeric/JacobianBasis.h b/Numeric/JacobianBasis.h
index 312808ab890eb5f65b6807f2fe7a1ba23b1974bd..dc9c8b598719daa63bc26e8b14eff43de96f92fd 100644
--- a/Numeric/JacobianBasis.h
+++ b/Numeric/JacobianBasis.h
@@ -16,30 +16,17 @@ class GradientBasis {
public:
fullMatrix<double> gradShapeMatX, gradShapeMatY, gradShapeMatZ;
- private:
- const bezierBasis *_bezier;
-
public :
GradientBasis(int tag, int order);
int getNumSamplingPoints() const {return gradShapeMatX.size1();}
int getNumMapNodes() const {return gradShapeMatX.size2();}
- const bezierBasis* getBezierBasis() const {return _bezier;}
void getGradientsFromNodes(const fullMatrix<double> &nodes,
fullMatrix<double> *dxyzdX,
fullMatrix<double> *dxyzdY,
fullMatrix<double> *dxyzdZ) const;
- inline double getgradShapeMatX(int i, int j) const {return gradShapeMatX(i, j);}
- inline double getgradShapeMatY(int i, int j) const {return gradShapeMatY(i, j);}
- inline double getgradShapeMatZ(int i, int j) const {return gradShapeMatZ(i, j);}
- inline void lag2Bez(const fullVector<double> &jac, fullVector<double> &bez) const {
- _bezier->matrixLag2Bez.mult(jac,bez);
- }
- inline void lag2Bez(const fullMatrix<double> &jac, fullMatrix<double> &bez) const {
- _bezier->matrixLag2Bez.mult(jac,bez);
- }
};
diff --git a/Numeric/MetricBasis.cpp b/Numeric/MetricBasis.cpp
index 865f31aee3e7d71290fd858761cd7f042c6228d4..a77e4e8c30b480396f414d73a636cf697d024191 100644
--- a/Numeric/MetricBasis.cpp
+++ b/Numeric/MetricBasis.cpp
@@ -12,9 +12,11 @@
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);
_gradients = BasisFactory::getGradientBasis(tag, metricOrder);
+ _bezier = BasisFactory::getBezierBasis(type, metricOrder);
int nSampPnts = _gradients->getNumSamplingPoints();
int nMapping = _gradients->getNumMapNodes();
@@ -143,7 +145,7 @@ void MetricCoefficient::interpolate(const double *uvw, double *minmaxQ)
return;
}
- int order = _gradients->getBezierBasis()->getOrder();
+ int order = _bezier->getOrder();
int dimSimplex;
fullMatrix<double> exponents;
@@ -277,6 +279,69 @@ void MetricCoefficient::interpolate(const double *uvw, double *minmaxQ)
}
}
+double MetricCoefficient::getBoundRmin(double tol)
+{
+ 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);
+ }
+ }
+
+ // 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);
+ }
+ else {
+
+ }*/
+
+
+ /*compute RminBez
+ if RminLag-RminBez < tol : return RminBez
+ return subdivide*/
+
+}
int MetricCoefficient::metricOrder(int tag)
{
@@ -355,5 +420,5 @@ void MetricCoefficient::_computeBezCoeff()
{
_coefficientsBez.resize(_coefficientsLag.size1(),
_coefficientsLag.size2());
- _gradients->lag2Bez(_coefficientsLag, _coefficientsBez);
+ _bezier->matrixLag2Bez.mult(_coefficientsLag, _coefficientsBez);
}
diff --git a/Numeric/MetricBasis.h b/Numeric/MetricBasis.h
index d1cc8a98fc29d3f3fcdc825e3047ebf6feb8329d..dcf1ab97d919af34e74651c113f720e6627a9a2e 100644
--- a/Numeric/MetricBasis.h
+++ b/Numeric/MetricBasis.h
@@ -15,6 +15,7 @@ class MetricCoefficient {
MElement *_element;
const JacobianBasis *_jacobian;
const GradientBasis *_gradients;
+ const bezierBasis *_bezier;
fullMatrix<double> _coefficientsLag, _coefficientsBez;
public:
@@ -22,7 +23,7 @@ class MetricCoefficient {
void getCoefficients(fullMatrix<double>&, bool bezier = true);
void interpolate(const double *uvw, double *minmaxQ);
- void getBoundRmin(double tol);
+ double getBoundRmin(double tol);
static int metricOrder(int tag);