From 75cd51b9aea1819700ab5aa8da64e1eac44bf5e8 Mon Sep 17 00:00:00 2001
From: Christophe Geuzaine <cgeuzaine@ulg.ac.be>
Date: Thu, 31 Mar 2016 08:54:50 +0000
Subject: [PATCH] alternative addPlanarFace from Wendy

---
 Geo/GEdgeLoop.h       |   2 +-
 Geo/GModel.cpp        |   7 ++
 Geo/GModel.h          |   3 +-
 Geo/GModelFactory.cpp | 215 +++++++++++++++++++++++++-----------------
 Geo/GModelFactory.h   |   8 ++
 5 files changed, 144 insertions(+), 91 deletions(-)

diff --git a/Geo/GEdgeLoop.h b/Geo/GEdgeLoop.h
index 54982ceed2..84e6f149bc 100644
--- a/Geo/GEdgeLoop.h
+++ b/Geo/GEdgeLoop.h
@@ -22,7 +22,7 @@ struct GEdgeSigned
     return (_sign != 1) ? ge->getBeginVertex() : ge->getEndVertex();
   }
   void print() const;
-  int getSign(){ return _sign; }
+  int getSign() const { return _sign; }
 };
 
 class GEdgeLoop
diff --git a/Geo/GModel.cpp b/Geo/GModel.cpp
index 95f0fa35ce..7b12183962 100644
--- a/Geo/GModel.cpp
+++ b/Geo/GModel.cpp
@@ -2946,6 +2946,13 @@ GFace* GModel::addPlanarFace (std::vector<std::vector<GEdge *> > edges)
   return 0;
 }
 
+GFace* GModel::addPlanarFace (std::vector<std::vector<GEdgeSigned> > edges)
+{
+  if(_factory)
+    return _factory->addPlanarFace(this, edges);
+  return 0;
+}
+
 GRegion* GModel::addVolume (std::vector<std::vector<GFace *> > faces)
 {
   if(_factory)
diff --git a/Geo/GModel.h b/Geo/GModel.h
index 96722dde8d..253172e5df 100644
--- a/Geo/GModel.h
+++ b/Geo/GModel.h
@@ -121,7 +121,7 @@ class GModel
 
   // store the parent's pointer back into MSubElements (replacing numeric id)
   void _storeParentsInSubElements(std::map< int, std::vector<MElement* > >& map);
-  
+
   // Discrete Entities have to have their mesh moved to a geometry container
   void _createGeometryOfDiscreteEntities(bool force=false);
 
@@ -523,6 +523,7 @@ class GModel
   std::vector<GFace *> addRuledFaces(std::vector<std::vector<GEdge *> > edges);
   GFace *addFace(std::vector<GEdge *> edges, std::vector< std::vector<double > > points);
   GFace *addPlanarFace(std::vector<std::vector<GEdge *> > edges);
+  GFace *addPlanarFace (std::vector<std::vector<GEdgeSigned> > edges);
   GFace *add2Drect(double x0, double y0, double dx, double dy);
   GFace *add2Dellips(double xc, double yc, double rx, double ry);
 
diff --git a/Geo/GModelFactory.cpp b/Geo/GModelFactory.cpp
index 783d716d20..f738f488d4 100644
--- a/Geo/GModelFactory.cpp
+++ b/Geo/GModelFactory.cpp
@@ -66,99 +66,25 @@ GEdge *GeoFactory::addLine(GModel *gm, GVertex *start, GVertex *end)
 
 GFace *GeoFactory::addPlanarFace(GModel *gm, std::vector< std::vector<GEdge *> > edges)
 {
-
-  //create line loops
-  int nLoops = edges.size();
-  std::vector<EdgeLoop *> vecLoops;
-  for (int i=0; i< nLoops; i++){
-    int ne=(int)edges[i].size();
-    List_T *temp = List_Create(ne, ne, sizeof(int));
-    for(int j = 0; j < ne; j++){
-      GEdge *ge = edges[i][j];
-      int numEdge = ge->tag();
-      //create curve if it does not exist
-      Curve *c = FindCurve(numEdge);
-      if(!c){
-	GVertex *gvb = ge->getBeginVertex();
-	GVertex *gve = ge->getEndVertex();
-	Vertex *vertb = FindPoint(abs(gvb->tag()));
-	Vertex *verte = FindPoint(abs(gve->tag()));
-	if (!vertb){
-	  vertb = Create_Vertex(gvb->tag(), gvb->x(), gvb->y(), gvb->z(),
-				gvb->prescribedMeshSizeAtVertex(), 1.0);
-	  Tree_Add(gm->getGEOInternals()->Points, &vertb);
-	  vertb->Typ = MSH_POINT;
-	  vertb->Num = gvb->tag();
-	 }
-	if (!verte){
-	  verte = Create_Vertex(gve->tag(), gve->x(), gve->y(), gve->z(),
-				gve->prescribedMeshSizeAtVertex(), 1.0);
-	  Tree_Add(gm->getGEOInternals()->Points, &verte);
-	  verte->Typ = MSH_POINT;
-	  verte->Num = gve->tag();
-	}
-
-	if (ge->geomType() == GEntity::Line){
-	  c = Create_Curve(numEdge, MSH_SEGM_LINE, 1, NULL, NULL, -1, -1, 0., 1.);
-	}
-	else if (ge->geomType() == GEntity::DiscreteCurve){
-	  c = Create_Curve(numEdge, MSH_SEGM_DISCRETE, 1, NULL, NULL, -1, -1, 0., 1.);
-	}
-	else if(ge->geomType() == GEntity::CompoundCurve){
-	  c = Create_Curve(numEdge, MSH_SEGM_COMPOUND, 1, NULL, NULL, -1, -1, 0., 1.);
-	  std::vector<GEdge*> gec = ((GEdgeCompound*)ge)->getCompounds();
-	  for(unsigned int i = 0; i < gec.size(); i++)
-	    c->compound.push_back(gec[i]->tag());
-	}
-	else{
-	  c = Create_Curve(numEdge, MSH_SEGM_DISCRETE, 1, NULL, NULL, -1, -1, 0., 1.);
-	}
-
-	c->Control_Points = List_Create(2, 1, sizeof(Vertex *));
-	List_Add(c->Control_Points, &vertb);
-	List_Add(c->Control_Points, &verte);
-	c->beg = vertb;
-	c->end = verte;
-	End_Curve(c);
-
-	Tree_Add(gm->getGEOInternals()->Curves, &c);
-	CreateReversedCurve(c);
-      }
-      List_Add(temp, &numEdge);
-    }
-
-    int numl = gm->getMaxElementaryNumber(1) + i;
-    while (FindEdgeLoop(numl)){
-      numl++;
-      if (!FindEdgeLoop(numl)) break;
+  std::vector<std::vector<GEdgeSigned> > orientedEdges;
+  orientedEdges.reserve(edges.size());
+  for (size_t i=0; i< edges.size(); i++){
+    std::vector<GEdge *> loop = edges[i];
+    orientedEdges.push_back(std::vector<GEdgeSigned>());
+    std::vector<GEdgeSigned> &orientedLoop = orientedEdges.back();
+    orientedLoop.reserve(loop.size());
+    for (size_t j=0; j< loop.size(); j++){
+      GEdge *edge = loop[j];
+      orientedLoop.push_back(GEdgeSigned(1, edge));
     }
-    sortEdgesInLoop(numl, temp, true);
-    EdgeLoop *l = Create_EdgeLoop(numl, temp);
-    vecLoops.push_back(l);
-    Tree_Add(gm->getGEOInternals()->EdgeLoops, &l);
-    l->Num = numl;
-    List_Delete(temp);
   }
 
-  //create surface
-  int numf  = gm->getMaxElementaryNumber(2)+1;
-  Surface *s = Create_Surface(numf, MSH_SURF_PLAN);
-  List_T *temp = List_Create(nLoops, nLoops, sizeof(int));
-  for (int i = 0; i < nLoops; i++){
-    int numl = vecLoops[i]->Num;
-    List_Add(temp, &numl);
-  }
-
-  setSurfaceGeneratrices(s, temp);
-  List_Delete(temp);
-  End_Surface(s);
-  Tree_Add(gm->getGEOInternals()->Surfaces, &s);
-
-  //gmsh surface
-  GFace *gf = new gmshFace(gm,s);
-  gm->add(gf);
+  return _addPlanarFace(gm, orientedEdges, true);
+}
 
-  return gf;
+GFace *GeoFactory::addPlanarFace(GModel *gm, const std::vector<std::vector<GEdgeSigned> > &edges)
+{
+  return _addPlanarFace(gm, edges, false);
 }
 
 GRegion* GeoFactory::addVolume (GModel *gm, std::vector<std::vector<GFace *> > faces)
@@ -435,6 +361,105 @@ void GeoFactory::healGeometry(GModel *gm, double tolerance)
   GModel::setCurrent(current);
 }
 
+GFace *GeoFactory::_addPlanarFace(GModel *gm, const std::vector<std::vector<GEdgeSigned> > &edges, bool orientEdges)
+{
+
+  //create line loops
+  int nLoops = edges.size();
+  std::vector<EdgeLoop *> vecLoops;
+  for (int i=0; i< nLoops; i++){
+    int ne=(int)edges[i].size();
+    List_T *temp = List_Create(ne, ne, sizeof(int));
+    for(int j = 0; j < ne; j++){
+      const GEdgeSigned &signedEdge = edges[i][j];
+      GEdge *ge = signedEdge.ge;
+      int numEdge = ge->tag();
+      //create curve if it does not exist
+      Curve *c = FindCurve(numEdge);
+      if(!c){
+	GVertex *gvb = ge->getBeginVertex();
+	GVertex *gve = ge->getEndVertex();
+	Vertex *vertb = FindPoint(abs(gvb->tag()));
+	Vertex *verte = FindPoint(abs(gve->tag()));
+	if (!vertb){
+	  vertb = Create_Vertex(gvb->tag(), gvb->x(), gvb->y(), gvb->z(),
+				gvb->prescribedMeshSizeAtVertex(), 1.0);
+	  Tree_Add(gm->getGEOInternals()->Points, &vertb);
+	  vertb->Typ = MSH_POINT;
+	  vertb->Num = gvb->tag();
+	 }
+	if (!verte){
+	  verte = Create_Vertex(gve->tag(), gve->x(), gve->y(), gve->z(),
+				gve->prescribedMeshSizeAtVertex(), 1.0);
+	  Tree_Add(gm->getGEOInternals()->Points, &verte);
+	  verte->Typ = MSH_POINT;
+	  verte->Num = gve->tag();
+	}
+
+	if (ge->geomType() == GEntity::Line){
+	  c = Create_Curve(numEdge, MSH_SEGM_LINE, 1, NULL, NULL, -1, -1, 0., 1.);
+	}
+	else if (ge->geomType() == GEntity::DiscreteCurve){
+	  c = Create_Curve(numEdge, MSH_SEGM_DISCRETE, 1, NULL, NULL, -1, -1, 0., 1.);
+	}
+	else if(ge->geomType() == GEntity::CompoundCurve){
+	  c = Create_Curve(numEdge, MSH_SEGM_COMPOUND, 1, NULL, NULL, -1, -1, 0., 1.);
+	  std::vector<GEdge*> gec = ((GEdgeCompound*)ge)->getCompounds();
+	  for(unsigned int i = 0; i < gec.size(); i++)
+	    c->compound.push_back(gec[i]->tag());
+	}
+	else{
+	  c = Create_Curve(numEdge, MSH_SEGM_DISCRETE, 1, NULL, NULL, -1, -1, 0., 1.);
+	}
+
+	c->Control_Points = List_Create(2, 1, sizeof(Vertex *));
+	List_Add(c->Control_Points, &vertb);
+	List_Add(c->Control_Points, &verte);
+	c->beg = vertb;
+	c->end = verte;
+	End_Curve(c);
+
+	Tree_Add(gm->getGEOInternals()->Curves, &c);
+	CreateReversedCurve(c);
+      }
+      int signedNumEdge = numEdge*signedEdge.getSign();
+      List_Add(temp, &signedNumEdge);
+    }
+
+    int numl = gm->getMaxElementaryNumber(1) + i;
+    while (FindEdgeLoop(numl)){
+      numl++;
+      if (!FindEdgeLoop(numl)) break;
+    }
+    sortEdgesInLoop(numl, temp, orientEdges);
+    EdgeLoop *l = Create_EdgeLoop(numl, temp);
+    vecLoops.push_back(l);
+    Tree_Add(gm->getGEOInternals()->EdgeLoops, &l);
+    l->Num = numl;
+    List_Delete(temp);
+  }
+
+  //create surface
+  int numf  = gm->getMaxElementaryNumber(2)+1;
+  Surface *s = Create_Surface(numf, MSH_SURF_PLAN);
+  List_T *temp = List_Create(nLoops, nLoops, sizeof(int));
+  for (int i = 0; i < nLoops; i++){
+    int numl = vecLoops[i]->Num;
+    List_Add(temp, &numl);
+  }
+
+  setSurfaceGeneratrices(s, temp);
+  List_Delete(temp);
+  End_Surface(s);
+  Tree_Add(gm->getGEOInternals()->Surfaces, &s);
+
+  //gmsh surface
+  GFace *gf = new gmshFace(gm,s);
+  gm->add(gf);
+
+  return gf;
+}
+
 #if defined(HAVE_OCC)
 #include "OCCIncludes.h"
 #include "GModelIO_OCC.h"
@@ -1418,6 +1443,12 @@ GFace *OCCFactory::addPlanarFace(GModel *gm, std::vector< std::vector<GEdge *> >
   return gm->_occ_internals->addFaceToModel(gm, TopoDS::Face(aResult));
 }
 
+GFace *OCCFactory::addPlanarFace(GModel *gm, const std::vector<std::vector<GEdgeSigned> > &edges)
+{
+  Msg::Error("addPlanarFace with oriented edges not implemented yet for OCCFactory");
+  return 0;
+}
+
 GEntity *OCCFactory::addPipe(GModel *gm, GEntity *base, std::vector<GEdge *> wire)
 {
   BRepBuilderAPI_MakeWire wire_maker;
@@ -1497,6 +1528,12 @@ GFace* SGEOMFactory::addPlanarFace(GModel *gm, std::vector<std::vector<GEdge *>
   return 0;
 }
 
+GFace* SGEOMFactory::addPlanarFace(GModel *gm, const std::vector<std::vector<GEdgeSigned> > &edges)
+{
+  Msg::Error("addPlanarFace with oriented edges not implemented yet for SGEOMFactory");
+  return 0;
+}
+
 GRegion* SGEOMFactory::addVolume(GModel *gm, std::vector<std::vector<GFace *> > faces)
 {
   Msg::Error("addVolume not implemented yet for SGEOMFactory");
diff --git a/Geo/GModelFactory.h b/Geo/GModelFactory.h
index 1aee7a9b3b..00cff49b0a 100644
--- a/Geo/GModelFactory.h
+++ b/Geo/GModelFactory.h
@@ -13,6 +13,7 @@
 class GEntity;
 class GVertex;
 class GEdge;
+class GEdgeSigned;
 class GFace;
 class GRegion;
 class GModel;
@@ -31,6 +32,7 @@ class GModelFactory {
                              double lc) = 0;
   virtual GEdge *addLine(GModel *, GVertex *v1, GVertex *v2) = 0;
   virtual GFace *addPlanarFace(GModel *gm, std::vector<std::vector<GEdge *> > edges) = 0;
+  virtual GFace *addPlanarFace(GModel *gm, const std::vector<std::vector<GEdgeSigned> > &edges) = 0;
   virtual GRegion*addVolume(GModel *gm, std::vector<std::vector<GFace *> > faces) = 0;
   virtual GEdge *addCircleArc(GModel *gm, GVertex *start, GVertex *center, GVertex *end)
   {
@@ -219,12 +221,16 @@ class GeoFactory : public GModelFactory {
   GVertex *addVertex(GModel *gm,double x, double y, double z, double lc);
   GEdge *addLine(GModel *gm,GVertex *v1, GVertex *v2);
   GFace *addPlanarFace(GModel *gm, std::vector<std::vector<GEdge *> > edges);
+  GFace *addPlanarFace(GModel *gm, const std::vector<std::vector<GEdgeSigned> > &edges);
   GRegion *addVolume(GModel *gm, std::vector<std::vector<GFace *> > faces);
   GEdge *addCircleArc(GModel *gm,GVertex *begin, GVertex *center, GVertex *end);
   std::vector<GFace *> addRuledFaces(GModel *gm, std::vector<std::vector<GEdge *> > edges);
   std::vector<GEntity*> extrudeBoundaryLayer(GModel *gm, GEntity *e, int nbLayers,
                                              double hLayers, int dir, int view);
   void healGeometry(GModel *gm, double tolerance = -1.);
+
+private:
+  GFace *_addPlanarFace(GModel *gm, const std::vector<std::vector<GEdgeSigned> > &edges, bool orientEdges);
 };
 
 #if defined(HAVE_OCC)
@@ -258,6 +264,7 @@ class OCCFactory : public GModelFactory {
   GFace *addFace(GModel *gm, std::vector<GEdge *> edges,
                  std::vector< std::vector<double > > points);
   GFace *addPlanarFace(GModel *gm, std::vector<std::vector<GEdge *> > edges);
+  GFace *addPlanarFace(GModel *gm, const std::vector<std::vector<GEdgeSigned> > &edges);
   GFace *add2Drect(GModel *gm,double x0, double y0, double dx, double dy);
   GFace *add2Dellips(GModel *gm,double xc, double yc, double rx, double ry);
   GRegion *addVolume(GModel *gm, std::vector<std::vector<GFace *> > faces);
@@ -293,6 +300,7 @@ class SGEOMFactory : public GModelFactory {
   GVertex *addVertex(GModel *gm,double x, double y, double z, double lc);
   GEdge *addLine(GModel *gm,GVertex *v1, GVertex *v2);
   GFace *addPlanarFace(GModel *gm, std::vector<std::vector<GEdge *> > edges);
+  GFace *addPlanarFace(GModel *gm, const std::vector<std::vector<GEdgeSigned> > &edges);
   GRegion *addVolume(GModel *gm, std::vector<std::vector<GFace *> > faces);
   void healGeometry(GModel *gm, double tolerance = -1.);
 };
-- 
GitLab