Fixed orientation of surface meshes with boundary layers

Mesh/meshGFace.cpp

@@ -2814,87 +2814,127 @@ void partitionAndRemesh(GFaceCompound *gf)
 #endif
 }
 
-void orientMeshGFace::operator()(GFace *gf)
+static bool getGFaceNormalFromVert(GFace *gf, MElement *el, SVector3 &nf)
 {
-  if(!gf->getNumMeshElements()) return;
-  if(gf->geomType() == GEntity::ProjectionFace) return;
-
-  gf->model()->setCurrentMeshEntity(gf);
-
-  if(gf->geomType() == GEntity::DiscreteSurface ||
-     gf->geomType() == GEntity::BoundaryLayerSurface ||
-     gf->geomType() == GEntity::CompoundSurface){
-    // don't do anything
-  }
-  else{
-    // in old versions we checked the orientation by comparing the orientation
-    // of a line element on the boundary w.r.t. its connected surface
-    // element. This is probably better than what follows, but
-    // * it failed when the 3D Delaunay changes the 1D mesh (since we don't
-    //    recover it yet)
-    // * it failed with OpenCASCADE geometries, where surface orientions do not
-    //   seem to be consistent with the orientation of the bounding edges
-
-    bool done = false;
-    // first, try to find an element with one vertex categorized on the
-    // surface and for which we have valid surface parametric
-    // coordinates
-    for(unsigned int i = 0; i < gf->getNumMeshElements(); i++){
-      MElement *e = gf->getMeshElement(i);
-      for(int j = 0; j < e->getNumVertices(); j++){
-        MVertex *v = e->getVertex(j);
+  bool found = false;
+  for (int iElV = 0; iElV < el->getNumVertices(); iElV++) {
+    MVertex *v = el->getVertex(iElV);
     SPoint2 param;
     if(v->onWhat() == gf && v->getParameter(0, param[0]) &&
        v->getParameter(1, param[1])) {
-          SVector3 nf = gf->normal(param);
-          SVector3 ne = e->getFace(0).normal();
-          if(dot(ne, nf) < 0){
-            Msg::Debug("Reversing orientation of mesh in face %d (param)", gf->tag());
-            for(unsigned int k = 0; k < gf->getNumMeshElements(); k++)
-              gf->getMeshElement(k)->reverse();
-          }
-          done = true;
+      nf = gf->normal(param);
+      found = true;
       break;
     }
   }
-      if(done) break;
+  return found;
 }
 
-    if(!done){
-      // if we could not find such an element, just try to evaluate the
-      // normal at the barycenter of an element on the surface
-      for(unsigned int i = 0; i < gf->getNumMeshElements(); i++){
-        MElement *e = gf->getMeshElement(i);
+static bool getGFaceNormalFromBary(GFace *gf, MElement *el, SVector3 &nf)
+{
   SPoint2 param(0., 0.);
   bool ok = true;
-        for(int j = 0; j < e->getNumVertices(); j++){
+  for (int j = 0; j < el->getNumVertices(); j++) {
     SPoint2 p;
     // FIXME: use inexact reparam because some vertices might not be
     // exactly on the surface after the 3D Delaunay
-          ok = reparamMeshVertexOnFace(e->getVertex(j), gf, p, false);
+    ok = reparamMeshVertexOnFace(el->getVertex(j), gf, p, false);
     if (!ok) break;
     param += p;
   }
   if (ok) {
-          param *= 1. / e->getNumVertices();
-          SVector3 nf = gf->normal(param);
+    param *= 1. / el->getNumVertices();
+    nf = gf->normal(param);
+  }
+  return ok;
+}
+
+static void getGFaceOrientation(GFace *gf, BoundaryLayerColumns *blc,
+                                bool existBL, bool fromVert,
+                                int &orientNonBL, int &orientBL)
+{
+  for (unsigned int iEl = 0; iEl < gf->getNumMeshElements(); iEl++) {
+    MElement *e = gf->getMeshElement(iEl);
+    const bool isBLEl = existBL &&
+                        (blc->_toFirst.find(e) != blc->_toFirst.end());
+    SVector3 nf;
+    if ((!isBLEl && orientNonBL == 0) || (isBLEl && orientBL == 0)) {       // Check only if orientation of BL/non-BL el. not already known
+      const bool found = fromVert ? getGFaceNormalFromVert(gf, e, nf) :
+                                    getGFaceNormalFromBary(gf, e, nf);
+      if (found) {
         SVector3 ne = e->getFace(0).normal();
-          if(dot(ne, nf) < 0){
-            Msg::Debug("Reversing 2 orientation of mesh in face %d (cog)", gf->tag());
-            for(unsigned int k = 0; k < gf->getNumMeshElements(); k++)
-              gf->getMeshElement(k)->reverse();
+        const int orient = (dot(ne, nf) > 0.) ? 1 : -1;
+        if (isBLEl) orientBL = orient;
+        else orientNonBL = orient;
       }
-          done = true;
-          break;
     }
+    if ((orientNonBL != 0) && (orientBL != 0)) break;                       // Stop when orientation found for non-BL and BL el.
   }
 }
 
-    if(!done)
+void orientMeshGFace::operator()(GFace *gf)
+{
+  if(!gf->getNumMeshElements()) return;
+  if(gf->geomType() == GEntity::ProjectionFace) return;
+
+  gf->model()->setCurrentMeshEntity(gf);
+
+  if(gf->geomType() == GEntity::DiscreteSurface ||
+     gf->geomType() == GEntity::BoundaryLayerSurface ||
+     gf->geomType() == GEntity::CompoundSurface){
+    // don't do anything
+  }
+  else {
+    // In old versions we checked the orientation by comparing the orientation
+    // of a line element on the boundary w.r.t. its connected surface
+    // element. This is probably better than what follows, but
+    // * it failed when the 3D Delaunay changes the 1D mesh (since we don't
+    //    recover it yet)
+    // * it failed with OpenCASCADE geometries, where surface orientions do not
+    //   seem to be consistent with the orientation of the bounding edges
+
+    // Now: orient surface elements w.r.t. normal to geometric model.
+    // Assumes that originally, orientation is consistent among boundary layer
+    // (BL) elements, and orientation is consistent among non-BL elements, but
+    // BL and non-BL elements can be oriented differently
+
+    // Determine whether there is a boundary layer (BL)
+    BoundaryLayerColumns *blc = gf->getColumns();
+    const bool existBL = !blc->_toFirst.empty();
+
+    // Get orientation of BL and non-BL elements.
+    // First, try to get normal to GFace from vertices.
+    // If it fails, try to get normal to GFace from element barycenter
+    int orientNonBL = 0, orientBL = existBL ? 0 : 1;
+    getGFaceOrientation(gf, blc, existBL, true, orientNonBL, orientBL);
+    if ((orientNonBL == 0) || (orientBL == 0))
+      getGFaceOrientation(gf, blc, existBL, false, orientNonBL, orientBL);
+
+    // Exit if could not determine orientation of both non-BL el. and BL el.
+    if ((orientNonBL == 0) && (orientBL == 0)) {
       Msg::Warning("Could not orient mesh in face %d", gf->tag());
+      return;
+    }
+
+    // Reverse BL and non-BL elements if needed
+    if (existBL) {                                                          // If there is a BL, test BL/non-BL elements
+      if ((orientNonBL == -1) || (orientBL == -1))
+        for (unsigned int iEl = 0; iEl < gf->getNumMeshElements(); iEl++) {
+          MElement *e = gf->getMeshElement(iEl);
+          if (blc->_toFirst.find(e) == blc->_toFirst.end()) {               // If el. outside of BL...
+            if (orientNonBL == -1) e->reverse();                            // ... reverse if needed
+          }
+          else                                                              // If el. in BL
+            if (orientBL == -1) e->reverse();                               // ... reverse if needed
+        }
+    }
+    else                                                                    // If no BL, reverse all elements if needed
+      if (orientNonBL == -1)
+        for (unsigned int iEl = 0; iEl < gf->getNumMeshElements(); iEl++)
+          gf->getMeshElement(iEl)->reverse();
   }
 
-  // apply user-specified mesh orientation constraints
+  // Apply user-specified mesh orientation constraints
   if(gf->meshAttributes.reverseMesh)
     for(unsigned int k = 0; k < gf->getNumMeshElements(); k++)
       gf->getMeshElement(k)->reverse();