Merge branch 'master' of https://gitlab.onelab.info/gmsh/gmsh

Geo/OCCFace.cpp

@@ -256,7 +256,7 @@ GPoint OCCFace::closestPoint(const SPoint3 &qp,
   GeomAPI_ProjectPointOnSurf proj(pnt, occface, umin, umax, vmin, vmax);
 
   if(!proj.NbPoints()) {
-    Msg::Warning("OCC Project Point on Surface FAIL");
+    Msg::Debug("OCC Project Point on Surface FAIL");
     GPoint gp(0, 0);
     gp.setNoSuccess();
     return gp;

Mesh/meshGFace.cpp

@@ -131,7 +131,7 @@ bool pointInsideParametricDomain(std::vector<SPoint2> &bnd, SPoint2 &p,
 
 void trueBoundary(const char *iii, GFace *gf, std::vector<SPoint2> &bnd)
 {
-  //    FILE* view_t = Fopen(iii,"w");
+  ///FILE* view_t = Fopen(iii,"w");
   //      fprintf(view_t,"View \"True Boundary\"{\n");
   std::vector<GEdge *> edg = gf->edges();
 
@@ -2732,17 +2732,44 @@ void deMeshGFace::operator()(GFace *gf)
   gf->correspondingVertices.clear();
 }
 
-/*
-static double TRIANGLE_VALIDITY(MTriangle *t)
+
+static double TRIANGLE_VALIDITY(GFace *gf, MTriangle *t)
 {
+  SPoint2 p1,p2,p3;
+  reparamMeshVertexOnFace(t->getVertex(0),gf, p1);
+  reparamMeshVertexOnFace(t->getVertex(1),gf, p2);
+  reparamMeshVertexOnFace(t->getVertex(2),gf, p3);
+  SVector3 N1 = gf->normal(p1);
+  SVector3 N2 = gf->normal(p2);
+  SVector3 N3 = gf->normal(p3);
+  SVector3 N = N1 + N2 + N3;
+  N.normalize();
+  SVector3 d1(t->getVertex(1)->x() - t->getVertex(0)->x(),
+	      t->getVertex(1)->y() - t->getVertex(0)->y(),
+	      t->getVertex(1)->z() - t->getVertex(0)->z());
+  SVector3 d2(t->getVertex(2)->x() - t->getVertex(0)->x(),
+	      t->getVertex(2)->y() - t->getVertex(0)->y(),
+	      t->getVertex(2)->z() - t->getVertex(0)->z());
+
+  SVector3 c = crossprod (d1,d2);
+  
+  return dot (N,c);
+   
 }
-static bool isMeshValid (GFace *gf){
 
+static bool isMeshValid (GFace *gf){
+  int invalid = 0;
   for (size_t i=0;i<gf->triangles.size(); i++){
+    double v = TRIANGLE_VALIDITY(gf, gf->triangles[i]);
+    if (v < 0)invalid ++;
   }
+  if (invalid == 0 || invalid == gf->triangles.size())return true;
+  //  printf("%d %d %d\n",gf->tag(),gf->triangles.size(),invalid);
+  return false;
+  
   
 }
-*/
+
 
 // for debugging, change value from -1 to -100;
 int debugSurface = -1; //-100;
@@ -2837,6 +2864,15 @@ void meshGFace::operator()(GFace *gf, bool print)
   if(gf->getNumMeshElements() == 0) {
     Msg::Warning("Surface %d consists of no elements", gf->tag());
   }
+
+  if (algoDelaunay2D(gf) && !isMeshValid (gf)){
+    Msg::Warning ("Delaunay based mesher failed on surface %d --> moving to meshadapt",gf->tag());
+    deMeshGFace killer;
+    killer (gf);
+    gf->setMeshingAlgo (1);
+    (*this)(gf,print);
+  }
+
 }
 
 static bool getGFaceNormalFromVert(GFace *gf, MElement *el, SVector3 &nf)