- better surf orientation, but still not great (sigh)

- boundary of combined entities in parser

CMakeLists.txt

@@ -66,7 +66,8 @@ set(GMSH_API
   ${CMAKE_CURRENT_BINARY_DIR}/Common/GmshConfig.h 
   ${CMAKE_CURRENT_BINARY_DIR}/Common/GmshVersion.h
   Common/Gmsh.h Common/GmshDefines.h Common/GmshMessage.h Common/VertexArray.h
-  Numeric/Numeric.h Numeric/Gauss.h Numeric/polynomialBasis.h Numeric/fullMatrix.h
+  Numeric/Numeric.h Numeric/Gauss.h Numeric/polynomialBasis.h
+    Numeric/JacobianBasis.h Numeric/fullMatrix.h
     Numeric/simpleFunction.h Numeric/cartesian.h
   Geo/GModel.h Geo/GEntity.h Geo/GPoint.h Geo/GVertex.h Geo/GEdge.h 
     Geo/GFace.h Geo/GRegion.h Geo/GEdgeLoop.h Geo/GEdgeCompound.h 

Geo/Geo.cpp

@@ -1872,7 +1872,16 @@ void SymmetryShapes(double A, double B, double C, double D, List_T *shapes)
     ReplaceAllDuplicates();
 }
 
-void BoundaryShapes(List_T *shapes, List_T *shapesBoundary)
+class ShapeLessThan{
+ public:
+  bool operator()(Shape *v1, Shape *v2) const
+  {
+    if(std::abs(v1->Num) < std::abs(v2->Num)) return true;
+    return false;
+  }
+};
+
+void BoundaryShapes(List_T *shapes, List_T *shapesBoundary, bool combined)
 {
   for(int i = 0; i < List_Nbr(shapes); i++) {
     Shape O;
@@ -1955,6 +1964,26 @@ void BoundaryShapes(List_T *shapes, List_T *shapesBoundary)
       break;
     }
   }
+
+  if(combined){
+    // compute boundary of the combined shapes
+    std::set<Shape*, ShapeLessThan> combined;
+    for(int i = 0; i < List_Nbr(shapesBoundary); i++){
+      Shape *s = (Shape*)List_Pointer(shapesBoundary, i);
+      std::set<Shape*, ShapeLessThan>::iterator it = combined.find(s);
+      if(it == combined.end())
+        combined.insert(s);
+      else
+        combined.erase(it);
+    }
+    List_T *tmp = List_Create(combined.size(), 10, sizeof(Shape));
+    for(std::set<Shape*, ShapeLessThan>::iterator it = combined.begin(); 
+        it != combined.end(); it++)
+      List_Add(tmp, *it);
+    List_Reset(shapesBoundary);
+    List_Copy(tmp, shapesBoundary);
+    List_Delete(tmp);
+  }
 }
 
 // Extrusion routines
@@ -2829,6 +2858,8 @@ static void ReplaceDuplicatePoints()
   }
   List_Delete(All);
 
+  // TODO: replace old points in physical groups
+
   Tree_Action(points2delete, Free_Vertex);
   Tree_Delete(points2delete);
   Tree_Delete(allNonDuplicatedPoints);
@@ -2934,6 +2965,8 @@ static void ReplaceDuplicateCurves()
   }
   List_Delete(All);
 
+  // TODO: replace old curves in physical groups
+
   Tree_Action(curves2delete, Free_Curve);
   Tree_Delete(curves2delete);
   Tree_Delete(allNonDuplicatedCurves);
@@ -3024,6 +3057,8 @@ static void ReplaceDuplicateSurfaces()
   }
   List_Delete(All);
 
+  // TODO: replace old surfaces in physical groups
+
   Tree_Action(surfaces2delete, Free_Surface);
   Tree_Delete(surfaces2delete);
   Tree_Delete(allNonDuplicatedSurfaces);

Geo/Geo.h

@@ -278,7 +278,7 @@ void DilatShapes(double X,double Y,double Z, double A, List_T *shapes);
 void RotateShapes(double Ax,double Ay,double Az,
                   double Px,double Py, double Pz, double alpha, List_T *shapes);
 void SymmetryShapes(double A,double B,double C, double D, List_T *shapes);
-void BoundaryShapes(List_T *shapes, List_T *shapesBoundary);
+void BoundaryShapes(List_T *shapes, List_T *shapesBoundary, bool combined);
 void CopyShape(int Type, int Num, int *New);
 void DeleteShape(int Type, int Num);
 void ColorShape(int Type, int Num, unsigned int Color);

Geo/MVertex.cpp

@@ -328,7 +328,6 @@ bool reparamMeshEdgeOnFace(MVertex *v1, MVertex *v2, GFace *gf,
 
 bool reparamMeshVertexOnFace(const MVertex *v, const GFace *gf, SPoint2 &param)
 {
-
   if (gf->geomType() == GEntity::CompoundSurface &&
       v->onWhat()->dim() < 2){
     GFaceCompound *gfc = (GFaceCompound*) gf;

Mesh/meshGFace.cpp

@@ -1814,6 +1814,30 @@ void orientMeshGFace::operator()(GFace *gf)
   //   do not seem to be consistent with the orientation of the
   //   bounding edges
 
+  // 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("Reverting orientation of mesh in face %d", gf->tag());
+          for(unsigned int k = 0; k < gf->getNumMeshElements(); k++)
+            gf->getMeshElement(k)->revert();
+        }
+        return;
+      }
+    }
+  }
+  
+  // 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.);

Parser/Gmsh.y

@@ -1649,7 +1649,10 @@ Transform :
         }
       }
       else if(!strcmp($1, "Boundary")){
-        BoundaryShapes($3, $$);
+        BoundaryShapes($3, $$, false);
+      }
+      else if(!strcmp($1, "CombinedBoundary")){
+        BoundaryShapes($3, $$, true);
       }
       else{
         yymsg(0, "Unknown command on multiple shapes: '%s'", $1);