diff --git a/Mesh/HighOrder.cpp b/Mesh/HighOrder.cpp
index 97431605d842a21c809b4fc80db247c30362ca9e..692be877de7a472b8b734d60912541e894fac969 100644
--- a/Mesh/HighOrder.cpp
+++ b/Mesh/HighOrder.cpp
@@ -26,12 +26,6 @@
#define SQU(a) ((a)*(a))
-static bool mappingIsInvertible(MTetrahedron *e)
-{
- if (e->getPolynomialOrder() == 1) return true;
- if (e->getPolynomialOrder() == 1) return e->distoShapeMeasure() > 0;
-}
-
// The aim here is to build a polynomial representation that consist
// in polynomial segments of equal length
@@ -536,7 +530,7 @@ static void reorientQuadPoints(std::vector<MVertex*> &vtcs, int orientation,
else if (p1 == 0 && p2 == 3){
for (int i = 4+4*nbP-1; i >= 4+3*nbP; i--) tmp[index++] = vtcs[i];
}
- else printf("ouuls\n");
+ else Msg::Error("ouuls");
}
for (int i=0;i<index;i++)vtcs[start+4+i] = tmp[i];
start += (4+index);
@@ -856,8 +850,6 @@ static void setHighOrder(GRegion *gr, edgeContainer &edgeVertices,
int nPts = 1, highOrderSmoother *displ2D = 0,
highOrderSmoother *displ3D = 0)
{
- int nbCorr = 0;
-
std::vector<MTetrahedron*> tetrahedra2;
for(unsigned int i = 0; i < gr->tetrahedra.size(); i++){
MTetrahedron *t = gr->tetrahedra[i];
@@ -871,33 +863,18 @@ static void setHighOrder(GRegion *gr, edgeContainer &edgeVertices,
else{
getFaceVertices(gr, t, vf, faceVertices, edgeVertices, linear, nPts);
ve.insert(ve.end(), vf.begin(), vf.end());
- MTetrahedronN incpl(t->getVertex(0), t->getVertex(1), t->getVertex(2), t->getVertex(3),
- ve, nPts + 1);
+ MTetrahedronN incpl(t->getVertex(0), t->getVertex(1), t->getVertex(2),
+ t->getVertex(3), ve, nPts + 1);
getRegionVertices(gr, &incpl, t, vr, linear, nPts);
ve.insert(ve.end(), vr.begin(), vr.end());
MTetrahedron* n = new MTetrahedronN(t->getVertex(0), t->getVertex(1),
t->getVertex(2), t->getVertex(3), ve, nPts + 1);
- if (!mappingIsInvertible(n)){
- Msg::Warning("Found invalid curved volume element (# %d in list)", i);
- }
tetrahedra2.push_back(n);
}
delete t;
}
gr->tetrahedra = tetrahedra2;
- std::vector<int> invalid;
- if (nbCorr != 0) {
- for(unsigned int i = 0; i < gr->tetrahedra.size(); i++)
- if (!mappingIsInvertible(gr->tetrahedra[i])) invalid.push_back(i);
- if (invalid.size()) {
- Msg::Warning("We have %d invalid elements remaining", (int)invalid.size());
- std::vector<int>::iterator iIter = invalid.begin();
- for (; iIter != invalid.end(); ++iIter)
- Msg::Warning("%d", *iIter);
- }
- }
-
std::vector<MHexahedron*> hexahedra2;
for(unsigned int i = 0; i < gr->hexahedra.size(); i++){
MHexahedron *h = gr->hexahedra[i];
@@ -928,14 +905,14 @@ static void setHighOrder(GRegion *gr, edgeContainer &edgeVertices,
else {
getFaceVertices(gr, h, vf, faceVertices, edgeVertices, linear, nPts);
ve.insert(ve.end(), vf.begin(), vf.end());
- MHexahedronN incpl(h->getVertex(0), h->getVertex(1), h->getVertex(2), h->getVertex(3),
- h->getVertex(4), h->getVertex(5), h->getVertex(6), h->getVertex(7),
- ve, nPts + 1);
+ MHexahedronN incpl(h->getVertex(0), h->getVertex(1), h->getVertex(2),
+ h->getVertex(3), h->getVertex(4), h->getVertex(5),
+ h->getVertex(6), h->getVertex(7), ve, nPts + 1);
getRegionVertices(gr, &incpl, h, vr, linear, nPts);
ve.insert(ve.end(), vr.begin(), vr.end());
- MHexahedron* n = new MHexahedronN(h->getVertex(0), h->getVertex(1), h->getVertex(2), h->getVertex(3),
- h->getVertex(4), h->getVertex(5), h->getVertex(6), h->getVertex(7),
- ve, nPts + 1);
+ MHexahedron* n = new MHexahedronN(h->getVertex(0), h->getVertex(1), h->getVertex(2),
+ h->getVertex(3), h->getVertex(4), h->getVertex(5),
+ h->getVertex(6), h->getVertex(7), ve, nPts + 1);
hexahedra2.push_back(n);
}
delete h;
@@ -1226,8 +1203,8 @@ void SetOrderN(GModel *m, int order, bool linear, bool incomplete)
// we do that model face by model face
std::vector<MElement*> bad;
double worst;
- // printJacobians(m, "smoothness_b.pos");
+ // printJacobians(m, "smoothness_b.pos");
// m->writeMSH("RAW.msh");
if (displ2D){
@@ -1253,10 +1230,9 @@ void SetOrderN(GModel *m, int order, bool linear, bool incomplete)
if(displ2D) delete displ2D;
if(displ3D) delete displ3D;
- // printJacobians(m, "smoothness.pos");
-
double t2 = Cpu();
+ // printJacobians(m, "smoothness.pos");
// m->writeMSH("SMOOTHED.msh");
if (!linear){
diff --git a/Mesh/highOrderTools.cpp b/Mesh/highOrderTools.cpp
index ea884c357ece57e0a21819ff511f9ed732b1230c..f017ceeeffd70e6616d8837adf0bc4749fba41a0 100644
--- a/Mesh/highOrderTools.cpp
+++ b/Mesh/highOrderTools.cpp
@@ -160,7 +160,8 @@ double highOrderTools::applySmoothingTo(GFace *gf, double tres, bool mixed)
}
-void highOrderTools::ensureMinimumDistorsion (double threshold){
+void highOrderTools::ensureMinimumDistorsion (double threshold)
+{
std::vector<MElement*> v;
if (_dim == 2){
for (GModel::fiter fit = _gm->firstFace(); fit != _gm->lastFace(); ++fit) {
@@ -663,20 +664,18 @@ double highOrderTools::apply_incremental_displacement (double max_incr,
}
// uncurve elements that are invalid
-void highOrderTools::ensureMinimumDistorsion (std::vector<MElement*> &all,
- double threshold){
- int num = 0;
- while(1){
+void highOrderTools::ensureMinimumDistorsion(std::vector<MElement*> &all,
+ double threshold)
+{
+ for(int tries = 0; tries < 100; tries++){
double minD;
std::vector<MElement*> disto;
getDistordedElements(all, threshold, disto, minD);
- // if (num == disto.size())break;
- if (!disto.size())break;
- num = disto.size();
- Msg::Info("Fixing %d bad curved elements (worst disto %g)",disto.size(),minD);
- for (int i=0;i<disto.size();i++){
- ensureMinimumDistorsion(disto[i],threshold);
- }
+ if (disto.empty()) break;
+ Msg::Info("Fixing %d bad curved elements (worst disto %g)", disto.size(), minD);
+ for (int i = 0; i < disto.size(); i++){
+ ensureMinimumDistorsion(disto[i], threshold);
+ }
}
}