From d6f1c7d85fc6ce29332ba3f71bed2defd62da3e4 Mon Sep 17 00:00:00 2001
From: Emilie Marchandise <emilie.marchandise@uclouvain.be>
Date: Wed, 4 Apr 2012 10:22:06 +0000
Subject: [PATCH] iwork on extrudeBL

---
 Geo/GModel.cpp               | 77 ++++++++++++++++++++++++++++++++++--
 Geo/GModel.h                 |  3 +-
 Geo/GModelFactory.cpp        | 35 +++++++++++++++-
 Geo/Geo.cpp                  | 19 ++++++++-
 Mesh/CenterlineField.cpp     | 29 ++++++++------
 Mesh/meshGFaceExtruded.cpp   |  1 +
 Mesh/meshGRegionExtruded.cpp |  1 +
 7 files changed, 144 insertions(+), 21 deletions(-)

diff --git a/Geo/GModel.cpp b/Geo/GModel.cpp
index 9e066aad1d..428b00973f 100644
--- a/Geo/GModel.cpp
+++ b/Geo/GModel.cpp
@@ -2073,13 +2073,55 @@ void GModel::save(std::string fileName)
   CreateOutputFile(fileName, guess);
   GModel::setCurrent(temp);
 }
+GEdge* GModel::addCompoundEdge(std::vector<GEdge*> edges, int num){
 
-GFace* GModel::addCompoundFace(std::vector<GFace*> faces, int param, int typeS)
+  if (num ==-1) num =  getMaxElementaryNumber(1) + 1;
+  GEdgeCompound *gec = new GEdgeCompound(this, num, edges);
+  add(gec);
+
+  //create old geo format for the compound edge
+  //necessary for boundary layers
+  if(FindCurve(num)){
+    Msg::Error("Curve %d already exists", num);
+  }
+  else{
+    Curve *c = Create_Curve(num, MSH_SEGM_COMPOUND, 1, NULL, NULL, -1, -1, 0., 1.);
+    for(int i= 0; i< edges.size(); i++)
+      c->compound.push_back(edges[i]->tag());
+
+    // Curve *c = Create_Curve(num, MSH_SEGM_DISCRETE, 1,
+    // 			    NULL, NULL, -1, -1, 0., 1.);
+    // List_T *points = Tree2List(getGEOInternals()->Points);
+    // GVertex *gvb = gec->getBeginVertex();
+    // GVertex *gve = gec->getEndVertex();
+    // int nb = 2 ;
+    // c->Control_Points = List_Create(nb, 1, sizeof(Vertex *));
+    // for(int i = 0; i < List_Nbr(points); i++) {
+    //   Vertex *v;
+    //   List_Read(points, i, &v);
+    //   if (v->Num == gvb->tag()) {
+    // 	List_Add(c->Control_Points, &v);
+    // 	c->beg = v;
+    //   }
+    //   if (v->Num == gve->tag()) {
+    // 	List_Add(c->Control_Points, &v);
+    // 	c->end = v;
+    //   }
+    // }
+
+    End_Curve(c);
+    Tree_Add(getGEOInternals()->Curves, &c);
+    CreateReversedCurve(c);
+  }
+
+  return gec;
+}
+GFace* GModel::addCompoundFace(std::vector<GFace*> faces, int param, int split, int num)
 {
 #if defined(HAVE_SOLVER)
-  int num =  getMaxElementaryNumber(2) + 1;
+  if (num ==-1) num = getMaxElementaryNumber(2) + 1;
 
-  std::list<GFace*> comp(faces.begin(), faces.end());
+  std::list<GFace*> faces_comp(faces.begin(), faces.end());
   std::list<GEdge*> U0;
 
   GFaceCompound::typeOfCompound typ = GFaceCompound::HARMONIC_CIRCLE;
@@ -2091,7 +2133,34 @@ GFace* GModel::addCompoundFace(std::vector<GFace*> faces, int param, int typeS)
   if (param == 6) typ =  GFaceCompound::HARMONIC_SQUARE;
   if (param == 7) typ =  GFaceCompound::CONFORMAL_FE;
 
-  GFaceCompound *gfc = new GFaceCompound(this, num, comp, U0, typ, typeS);
+  GFaceCompound *gfc = new GFaceCompound(this, num, faces_comp, U0, typ, split);
+  
+  //create old geo format for the compound face
+  //necessary for boundary layers
+  if(FindSurface(num)){
+    Msg::Error("Surface %d already exists", num);
+  }
+  else{
+    Surface *s = Create_Surface(num, MSH_SURF_COMPOUND);
+    for(int i= 0; i< faces.size(); i++)
+      s->compound.push_back(faces[i]->tag());
+
+    // Surface *s = Create_Surface(num, MSH_SURF_DISCRETE);
+    // std::list<GEdge*> edges = gfc->edges();
+    // s->Generatrices = List_Create(edges.size(), 1, sizeof(Curve *));
+    // List_T *curves = Tree2List(_geo_internals->Curves);
+    // Curve *c;
+    // for(std::list<GEdge*>::iterator ite = edges.begin(); ite != edges.end(); ite++){
+    //   for(int i = 0; i < List_Nbr(curves); i++) {
+    // 	List_Read(curves, i, &c);
+    // 	if (c->Num == (*ite)->tag()) {
+    // 	  List_Add(s->Generatrices, &c);
+    // 	}
+    //   }
+    // }
+    
+    Tree_Add(_geo_internals->Surfaces, &s);
+  }
 
   add(gfc);
   return gfc;
diff --git a/Geo/GModel.h b/Geo/GModel.h
index f2812d1359..0c9b84f37d 100644
--- a/Geo/GModel.h
+++ b/Geo/GModel.h
@@ -456,7 +456,8 @@ class GModel
   void 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 *addCompoundFace(std::vector<GFace*> faces, int typeP, int typeS);
+  GEdge *addCompoundEdge(std::vector<GEdge*> edges, int num=-1);
+  GFace *addCompoundFace(std::vector<GFace*> faces, int type, int split, int num=-1);
   GRegion *addVolume(std::vector<std::vector<GFace*> > faces);
 
   // create solid geometry primitives using the factory
diff --git a/Geo/GModelFactory.cpp b/Geo/GModelFactory.cpp
index 765619e0e2..2f9e941d9a 100644
--- a/Geo/GModelFactory.cpp
+++ b/Geo/GModelFactory.cpp
@@ -297,8 +297,9 @@ GEntity* GeoFactory::extrudeBoundaryLayer(GModel *gm, GEntity *e, int nbLayers,
   ep.mesh.NbElmLayer.clear(); 
   ep.mesh.NbElmLayer.push_back(nbLayers);
   ep.mesh.ExtrudeMesh = true;
+  ep.geo.Source = e->tag();
 
-  int type  = BOUNDARY_LAYER; //TRANSLATE;
+  int type  = BOUNDARY_LAYER; 
   double T0=0., T1=0., T2=0.;
   double A0=0., A1=0., A2=0.;
   double X0=0., X1=0., X2=0.,alpha=0.;
@@ -321,8 +322,25 @@ GEntity* GeoFactory::extrudeBoundaryLayer(GModel *gm, GEntity *e, int nbLayers,
   else if (e->dim() == 2){
     ((GFace*)e)->meshAttributes.extrude = &ep;
     Surface *s = FindSurface(e->tag());
-    if(!s) printf("surface %d not found \n", e->tag());
+
+    Msg::Info("Geo internal model has:");
+    List_T *points = Tree2List(gm->getGEOInternals()->Points);
+    List_T *curves = Tree2List(gm->getGEOInternals()->Curves);
+    List_T *surfaces = Tree2List(gm->getGEOInternals()->Surfaces);
+    Msg::Info("%d Vertices", List_Nbr(points));
+    Msg::Info("%d Edges", List_Nbr(curves));
+    Msg::Info("%d Faces", List_Nbr(surfaces));
+    for(int i = 0; i < List_Nbr(surfaces); i++) {
+      Surface *s;
+      List_Read(surfaces, i, &s);
+      printf("surface %d \n", s->Num);
+    }
+    if(!s) {
+      printf("surface %d NOT found \n", e->tag());
+      exit(1);
+    }
     else printf("surface %d found type =%d\n", e->tag(), s->Typ);
+
     shape.Num = s->Num;
     shape.Type = s->Typ;
   }
@@ -338,6 +356,19 @@ GEntity* GeoFactory::extrudeBoundaryLayer(GModel *gm, GEntity *e, int nbLayers,
                 &ep,
                 list_out);
 
+  Msg::Info("AFTER EXTRUDE Geo internal model has:");
+  List_T *points = Tree2List(gm->getGEOInternals()->Points);
+  List_T *curves = Tree2List(gm->getGEOInternals()->Curves);
+  List_T *surfaces = Tree2List(gm->getGEOInternals()->Surfaces);
+  Msg::Info("%d Vertices", List_Nbr(points));
+  Msg::Info("%d Edges", List_Nbr(curves));
+  Msg::Info("%d Faces", List_Nbr(surfaces));
+  for(int i = 0; i < List_Nbr(surfaces); i++) {
+    Surface *s;
+    List_Read(surfaces, i, &s);
+    printf("surface %d \n", s->Num);
+  }
+
   //create GEntities 
   gm->importGEOInternals();
 
diff --git a/Geo/Geo.cpp b/Geo/Geo.cpp
index fbbd0d4239..6a8dbdd33e 100644
--- a/Geo/Geo.cpp
+++ b/Geo/Geo.cpp
@@ -2489,7 +2489,7 @@ int Extrude_ProtudeSurface(int type, int is,
     return 0;
 
   Msg::Debug("Extrude Surface %d", is);
-
+    
   chapeau = DuplicateSurface(ps, false);
   chapeau->Extrude = new ExtrudeParams(COPIED_ENTITY);
   chapeau->Extrude->fill(type, T0, T1, T2, A0, A1, A2, X0, X1, X2, alpha);
@@ -2512,7 +2512,7 @@ int Extrude_ProtudeSurface(int type, int is,
     c->Extrude->geo.Source = c2->Num;
     if(e)
       c->Extrude->mesh = e->mesh;
-  }
+  }   
 
   // FIXME: this is a really ugly hack for backward compatibility, so
   // that we don't screw up the old .geo files too much. (Before
@@ -2639,6 +2639,7 @@ int Extrude_ProtudeSurface(int type, int is,
   // this is done only for backward compatibility with the old
   // numbering scheme
   Tree_Suppress(GModel::current()->getGEOInternals()->Surfaces, &chapeau);
+
   chapeau->Num = NEWSURFACE();
   chapeau->meshMaster = chapeau->Num;
   GModel::current()->getGEOInternals()->MaxSurfaceNum = chapeau->Num;
@@ -2684,6 +2685,20 @@ void ExtrudeShapes(int type, List_T *list_in,
                    ExtrudeParams *e,
                    List_T *list_out)
 {
+
+   // Msg::Info("IN EXTRUDE SHAPE Geo internal model has:");
+   //  List_T *points = Tree2List(GModel::current()->getGEOInternals()->Points);
+   //  List_T *curves = Tree2List(GModel::current()->getGEOInternals()->Curves);
+   //  List_T *surfaces = Tree2List(GModel::current()->getGEOInternals()->Surfaces);
+   //  Msg::Info("%d Vertices", List_Nbr(points));
+   //  Msg::Info("%d Edges", List_Nbr(curves));
+   //  Msg::Info("%d Faces", List_Nbr(surfaces));
+   //  for(int i = 0; i < List_Nbr(surfaces); i++) {
+   //    Surface *s;
+   //    List_Read(surfaces, i, &s);
+   //    printf("surface %d \n", s->Num);
+   //  }
+
  
   for(int i = 0; i < List_Nbr(list_in); i++){
     Shape shape;
diff --git a/Mesh/CenterlineField.cpp b/Mesh/CenterlineField.cpp
index 9d1a9eab90..f936b76b4f 100644
--- a/Mesh/CenterlineField.cpp
+++ b/Mesh/CenterlineField.cpp
@@ -622,7 +622,7 @@ void Centerline::createSplitCompounds(){
   NR = current->getMaxElementaryNumber(3);
 
   // Remesh new faces (Compound Lines and Compound Surfaces)
-  Msg::Info("*** Starting parametrize compounds:");
+  Msg::Info("Centerline: create split compounds:");
 
   //Parametrize Compound Lines
   for (int i=0; i < NE; i++){
@@ -630,30 +630,34 @@ void Centerline::createSplitCompounds(){
     GEdge *pe = current->getEdgeByTag(i+1);//current edge
     e_compound.push_back(pe);
     int num_gec = NE+i+1;
-    Msg::Info("Parametrize Compound Line (%d) = %d discrete edge",
+    Msg::Info("Create Compound Line (%d) = %d discrete edge",
               num_gec, pe->tag());
-    GEdgeCompound *gec = new GEdgeCompound(current, num_gec, e_compound);
-    current->add(gec);
+    GEdge *gec = current->addCompoundEdge(e_compound,num_gec);
+    //GEdgeCompound *gec = new GEdgeCompound(current, num_gec, e_compound);
+    //current->add(gec);
   }
 
   // Parametrize Compound surfaces
   std::list<GEdge*> U0;
   for (int i=0; i < NF; i++){
-    std::list<GFace*> f_compound;
+    std::vector<GFace*> f_compound;
     GFace *pf =  current->getFaceByTag(i+1);//current face
     f_compound.push_back(pf);
     int num_gfc = NF+i+1;
-    Msg::Info("Parametrize Compound Surface (%d) = %d discrete face",
+    Msg::Info("Create Compound Surface (%d) = %d discrete face",
               num_gfc, pf->tag());
-    //GFaceCompound::typeOfCompound typ = GFaceCompound::CONVEX_CIRCLE; 
-    //GFaceCompound::typeOfCompound typ = GFaceCompound::HARMONIC_PLANE; 
-    GFaceCompound::typeOfCompound typ = GFaceCompound::CONFORMAL_SPECTRAL; 
-    GFaceCompound *gfc = new GFaceCompound(current, num_gfc, f_compound, U0,
-					   typ, 0);
+
+    GFace *gfc = current->addCompoundFace(f_compound, 1, 0, num_gfc); //1=conf_spectral 4=convex_circle
+    ////GFaceCompound::typeOfCompound typ = GFaceCompound::CONVEX_CIRCLE; 
+    ////GFaceCompound::typeOfCompound typ = GFaceCompound::HARMONIC_PLANE; 
+    //GFaceCompound::typeOfCompound typ = GFaceCompound::CONFORMAL_SPECTRAL; 
+    //GFaceCompound *gfc = new GFaceCompound(current, num_gfc, f_compound, U0,typ, 0);
+    //current->add(gfc);
+
     gfc->meshAttributes.recombine = recombine;
     gfc->addPhysicalEntity(100);
     current->setPhysicalName("newsurf", 2, 100);
-    current->add(gfc);
+
   }
 
 }
@@ -833,6 +837,7 @@ void Centerline::createClosedVolume(){
 void Centerline::extrudeBoundaryLayerWall(){
   
   printf("extrude boundary layer wall %d %g \n", nbElemLayer,  hLayer);
+
   for (int i= 0; i< NF; i++){
     GFace *gf = current->getFaceByTag(NF+i+1);//at this moment compound is not meshed yet exist yet
     current->setFactory("Gmsh");
diff --git a/Mesh/meshGFaceExtruded.cpp b/Mesh/meshGFaceExtruded.cpp
index 482801d3f7..62085422e3 100644
--- a/Mesh/meshGFaceExtruded.cpp
+++ b/Mesh/meshGFaceExtruded.cpp
@@ -262,6 +262,7 @@ int MeshExtrudedSurface(GFace *gf,
     GFace *from = gf->model()->getFaceByTag(std::abs(ep->geo.Source));
     if(!from){ 
       Msg::Error("Unknown source surface %d for extrusion", ep->geo.Source);
+      exit(1);
       return 0;
     }
     else if(from->geomType() != GEntity::DiscreteSurface &&
diff --git a/Mesh/meshGRegionExtruded.cpp b/Mesh/meshGRegionExtruded.cpp
index 060965577a..ebc8875d66 100644
--- a/Mesh/meshGRegionExtruded.cpp
+++ b/Mesh/meshGRegionExtruded.cpp
@@ -243,6 +243,7 @@ void meshGRegionExtruded::operator() (GRegion *gr)
   GFace *from = gr->model()->getFaceByTag(std::abs(ep->geo.Source));
   if(!from){
     Msg::Error("Unknown source surface %d for extrusion", ep->geo.Source);
+    exit(1);
     return;
   }
 
-- 
GitLab