diff --git a/Mesh/meshGFace.cpp b/Mesh/meshGFace.cpp
index 95248a9d9cf01c2851230b07b9655ca10bc5a612..1031141c8ae4f5d5b5969857485e1dd9651e1fc5 100644
--- a/Mesh/meshGFace.cpp
+++ b/Mesh/meshGFace.cpp
@@ -2814,6 +2814,64 @@ void partitionAndRemesh(GFaceCompound *gf)
#endif
}
+static bool getGFaceNormalFromVert(GFace *gf, MElement *el, SVector3 &nf)
+{
+ 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])) {
+ nf = gf->normal(param);
+ found = true;
+ break;
+ }
+ }
+ return found;
+}
+
+static bool getGFaceNormalFromBary(GFace *gf, MElement *el, SVector3 &nf)
+{
+ SPoint2 param(0., 0.);
+ bool ok = true;
+ 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(el->getVertex(j), gf, p, false);
+ if (!ok) break;
+ param += p;
+ }
+ if (ok) {
+ 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();
+ const int orient = (dot(ne, nf) > 0.) ? 1 : -1;
+ if (isBLEl) orientBL = orient;
+ else orientNonBL = orient;
+ }
+ }
+ if ((orientNonBL != 0) && (orientBL != 0)) break; // Stop when orientation found for non-BL and BL el.
+ }
+}
+
void orientMeshGFace::operator()(GFace *gf)
{
if(!gf->getNumMeshElements()) return;
@@ -2826,8 +2884,8 @@ void orientMeshGFace::operator()(GFace *gf)
gf->geomType() == GEntity::CompoundSurface){
// don't do anything
}
- else{
- // in old versions we checked the orientation by comparing the orientation
+ 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
@@ -2835,66 +2893,48 @@ void orientMeshGFace::operator()(GFace *gf)
// * 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);
- 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;
- break;
- }
- }
- if(done) break;
- }
-
- 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);
- SPoint2 param(0., 0.);
- bool ok = true;
- for(int j = 0; j < e->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);
- if(!ok) break;
- param += p;
- }
- if(ok){
- param *= 1. / e->getNumVertices();
- SVector3 nf = gf->normal(param);
- 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();
+ // 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
}
- done = true;
- break;
+ else // If el. in BL
+ if (orientBL == -1) e->reverse(); // ... reverse if needed
}
- }
}
-
- if(!done)
- Msg::Warning("Could not orient mesh in face %d", gf->tag());
+ 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();