diff --git a/Geo/discreteDiskFace.cpp b/Geo/discreteDiskFace.cpp
index b06f5660f59eab26f047973f4d493cd7b13aebc1..da0623e3f5914c39123db84ea8620e78d269b635 100644
--- a/Geo/discreteDiskFace.cpp
+++ b/Geo/discreteDiskFace.cpp
@@ -213,9 +213,9 @@ static bool orderVertices(const double &tot_length, const std::vector<MVertex*>
}
/*BUILDER*/
-discreteDiskFace::discreteDiskFace(int id, GFace *gf, triangulation* diskTriangulation,
+discreteDiskFace::discreteDiskFace(GFace *gf, triangulation* diskTriangulation,
int p, std::vector<GFace*> *CAD) :
- GFace(gf->model(),id), _parent (gf), _ddft(NULL), oct(NULL)
+ GFace(gf->model(),diskTriangulation->idNum), _parent (gf), _ddft(NULL), oct(NULL)
{
initialTriangulation = diskTriangulation;
std::vector<MElement*> mesh = diskTriangulation->tri;
@@ -275,7 +275,7 @@ discreteDiskFace::discreteDiskFace(int id, GFace *gf, triangulation* diskTriangu
discrete_triangles[i] = new MTriangle6 (vs);
}// end loop over triangles
- geoTriangulation = new triangulation(id,discrete_triangles,gf);
+ geoTriangulation = new triangulation(tag(),discrete_triangles,gf);
geoTriangulation->fillingHoles = diskTriangulation->fillingHoles;
allNodes = geoTriangulation->vert;
_totLength = geoTriangulation->bord.rbegin()->first;
@@ -288,17 +288,13 @@ discreteDiskFace::discreteDiskFace(int id, GFace *gf, triangulation* diskTriangu
if (!checkOrientationUV()){
Msg::Info("discreteDiskFace:: parametrization is not one-to-one; fixing "
"the discrete system.");
-
- //parametrize(true);
- //buildOct(CAD);
- optimize();
- buildOct(CAD);
-
- }
+ parametrize(true);
+ buildOct(CAD);
+ }
putOnView(true,false);
printParamMesh();
if(!checkOrientationUV()) Msg::Fatal("discreteDiskFace:: failing to fix the discrete system");
- else Msg::Info("Parameterization done :-)");
+ // else Msg::Info("Parameterization done :-)");
}
/*end BUILDER*/
@@ -320,10 +316,10 @@ void discreteDiskFace::buildOct(std::vector<GFace*> *CAD) const
oct = Octree_Create(maxElePerBucket, origin, ssize, discreteDiskFaceBB,
discreteDiskFaceCentroid, discreteDiskFaceInEle);
- _ddft = new discreteDiskFaceTriangle[discrete_triangles.size()-geoTriangulation->fillingHoles.size()];
+ _ddft = new discreteDiskFaceTriangle[discrete_triangles.size()/*-geoTriangulation->fillingHoles.size()*/];
int c = 0;
for(unsigned int i = 0; i < discrete_triangles.size(); ++i){
- if(geoTriangulation->fillingHoles.find(i)==geoTriangulation->fillingHoles.end()){
+ // if(geoTriangulation->fillingHoles.find(i)==geoTriangulation->fillingHoles.end()){
MElement *t = discrete_triangles[i];
discreteDiskFaceTriangle* my_ddft = &_ddft[c++];
my_ddft->p.resize(_N);
@@ -334,7 +330,7 @@ void discreteDiskFace::buildOct(std::vector<GFace*> *CAD) const
my_ddft->tri = t;
Octree_Insert(my_ddft, oct);
- }
+ // }
}
Octree_Arrange(oct);
}
@@ -466,19 +462,20 @@ bool discreteDiskFace::checkOrientationUV()
double p1[2] = {ct->p[0].x(), ct->p[0].y()};
double p2[2] = {ct->p[1].x(), ct->p[1].y()};
double p3[2] = {ct->p[2].x(), ct->p[2].y()};
- initial = robustPredicates::orient2d(p1, p2, p3);
- unsigned int i=1;
- for (; i<discrete_triangles.size()-geoTriangulation->fillingHoles.size(); i++){
+ unsigned int nbP = 0;
+ unsigned int nbM = 0;
+ for (unsigned int i=0; i<discrete_triangles.size(); i++){
ct = &_ddft[i];
p1[0] = ct->p[0].x(); p1[1] = ct->p[0].y();
p2[0] = ct->p[1].x(); p2[1] = ct->p[1].y();
p3[0] = ct->p[2].x(); p3[1] = ct->p[2].y();
current = robustPredicates::orient2d(p1, p2, p3);
- if(initial*current < 0.) {
- Msg::Error("Map %d of the atlas : Triangle UV %d has not the correct orientation (area %22.15E)",tag(),i+1,current);
- return false;
- break;
- }
+ if(current < 0.) nbM++;
+ else nbP++;
+ }
+ if (nbP*nbM){
+ Msg::Error("Map %d of the atlas : Triangles have different orientations (%d + / %d -)",tag(),nbP,nbM);
+ return false;
}
return true;
}
diff --git a/Geo/discreteDiskFace.h b/Geo/discreteDiskFace.h
index aa4908def9a7f21d9966566c2a754b5f9ca4a6f2..43f7621cb4ede103831dde56a625254119cf147f 100644
--- a/Geo/discreteDiskFace.h
+++ b/Geo/discreteDiskFace.h
@@ -80,7 +80,9 @@ class triangulation {
double aspectRatio()
{
double L = bord.rbegin()->first;
- return L / maxD;
+ if (L == 0.0)return 1.e22;
+ // printf("%12.5E %12.5E\n",L,maxD);
+ return maxD / L;
}
int genus()
@@ -222,7 +224,7 @@ class discreteDiskFace : public GFace {
void printParamMesh();
public:
- discreteDiskFace(int id, GFace *parent, triangulation* diskTriangulation,
+ discreteDiskFace(GFace *parent, triangulation* diskTriangulation,
int p=1, std::vector<GFace*> *CAD = NULL);
virtual ~discreteDiskFace();
void getTriangleUV(const double u,const double v,discreteDiskFaceTriangle **mt,
diff --git a/Geo/discreteFace.cpp b/Geo/discreteFace.cpp
index 1e5f01ae02b2ebc2cae5be547d051bbda510f939..f0c689d23f39bbfceca29444e0ca5b93ae3615b0 100644
--- a/Geo/discreteFace.cpp
+++ b/Geo/discreteFace.cpp
@@ -111,10 +111,11 @@ void discreteFace::secondDer(const SPoint2 ¶m,
void discreteFace::createGeometry()
{
checkAndFixOrientation();
- int order = 2;
+ int order = 1;
int nPart = 2;
- double eta = 2*3.14/7;
+ double eta = 5/(2.*3.14);
+
#if defined(HAVE_SOLVER) && defined(HAVE_ANN)
if (!_atlas.empty())return;
@@ -126,22 +127,25 @@ void discreteFace::createGeometry()
triangulation* init = new triangulation(-1, tem,this);
toSplit.push(init);
- if((toSplit.top())->genus()!=0 || (toSplit.top())->aspectRatio() < eta || (toSplit.top())->seamPoint){
+ if((toSplit.top())->genus()!=0 || (toSplit.top())->aspectRatio() > eta || (toSplit.top())->seamPoint){
+
while( !toSplit.empty()){
std::vector<triangulation*> part;
triangulation* tosplit = toSplit.top();
toSplit.pop();
-
+ // printf("genus %d ar %12.5E\n",tosplit->genus(), tosplit->aspectRatio());
+ // getchar();
split(tosplit,part,nPart);
delete tosplit;
for(unsigned int i=0; i<part.size(); i++){
- if(part[i]->genus()!=0 || part[i]->aspectRatio() < eta || part[i]->seamPoint)
+ if(part[i]->genus()!=0 || part[i]->aspectRatio() > eta || part[i]->seamPoint)
toSplit.push(part[i]);
else{
toParam.push_back(part[i]);
part[i]->idNum=id++;
+ // printf("part %d is OK\n", part[i]->idNum);
}
}// end for i
}// !.empty()
@@ -153,17 +157,17 @@ void discreteFace::createGeometry()
updateTopology(toParam);
for(unsigned int i=0; i<toParam.size(); i++){
- //printf("MAP(%d) : aspect ratio = %12.5E\n",i,toParam[i]->aspectRatio());
- //char name[256];
- //sprintf(name,"map%d.pos",i);
- //toParam[i]->print(name,i);
+ printf("MAP(%d) : aspect ratio = %12.5E\n",toParam[i]->idNum,toParam[i]->aspectRatio());
+ char name[256];
+ sprintf(name,"map%d.pos",i);
+ toParam[i]->print(name,i);
fillHoles(toParam[i]);
- //sprintf(name,"mapFilled%d.pos",i);
- //toParam[i]->print(name,i);
+ sprintf(name,"mapFilled%d.pos",i);
+ toParam[i]->print(name, toParam[i]->idNum);
}
for(unsigned int i=0; i<toParam.size(); i++){
std::vector<MElement*> mytri = toParam[i]->tri;
- discreteDiskFace *df = new discreteDiskFace (i,this,toParam[i], order,(_CAD.empty() ? NULL : &_CAD));
+ discreteDiskFace *df = new discreteDiskFace (this,toParam[i], order,(_CAD.empty() ? NULL : &_CAD));
df->replaceEdges(toParam[i]->my_GEdges);
_atlas.push_back(df);
}
diff --git a/contrib/HighOrderMeshOptimizer/OptHomMesh.cpp b/contrib/HighOrderMeshOptimizer/OptHomMesh.cpp
index 6c8534f6ebacaca6c784587f9dc97bddc402cc4c..afaa2d9dd4f87bd5055757dc42121a35af68541a 100644
--- a/contrib/HighOrderMeshOptimizer/OptHomMesh.cpp
+++ b/contrib/HighOrderMeshOptimizer/OptHomMesh.cpp
@@ -151,7 +151,7 @@ void Mesh::calcScaledNormalEl2D(const std::map<MElement*,GEntity*> &element2enti
if (scal < 0.) factor = -factor;
}
elNorm.scale(factor); // Re-scaling normal here is faster than an extra scaling operation on the Jacobian
-
+ // elNorm.setAll(1);
}
int Mesh::getFreeVertexStartIndex(MVertex* vert)