diff --git a/Geo/GModel.cpp b/Geo/GModel.cpp
index 428b00973f6eb4f93c96762503a672fa5d7f4abb..88b7ff256dcf99a416b515b8abbe67236475352b 100644
--- a/Geo/GModel.cpp
+++ b/Geo/GModel.cpp
@@ -2091,23 +2091,23 @@ GEdge* GModel::addCompoundEdge(std::vector<GEdge*> edges, int num){
// 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;
- // }
- // }
+ 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);
@@ -2145,19 +2145,34 @@ GFace* GModel::addCompoundFace(std::vector<GFace*> faces, int param, int split,
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);
- // }
- // }
- // }
+ std::list<GEdge*> edges = gfc->edges();
+
+ //Replace edges by their compounds
+ std::set<GEdge*> mySet;
+ std::list<GEdge*>::iterator it = edges.begin();
+ while(it != edges.end()){
+ if((*it)->getCompound()){
+ mySet.insert((*it)->getCompound());
+ }
+ else{
+ mySet.insert(*it);
+ }
+ ++it;
+ }
+ edges.clear();
+ edges.insert(edges.begin(), mySet.begin(), mySet.end());
+
+ 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);
}
diff --git a/Geo/GModelFactory.cpp b/Geo/GModelFactory.cpp
index da7c615ab6a5a944ccb6f437a34f80083283d652..9a7a5768db0cc871a862dcbe50031c169f232069 100644
--- a/Geo/GModelFactory.cpp
+++ b/Geo/GModelFactory.cpp
@@ -288,17 +288,17 @@ std::vector<GFace *> GeoFactory::addRuledFaces(GModel *gm,
GEntity* GeoFactory::extrudeBoundaryLayer(GModel *gm, GEntity *e, int nbLayers, double hLayer, int dir, int view)
{
- ExtrudeParams ep;
- ep.mesh.BoundaryLayerIndex = dir;
- ep.mesh.ViewIndex = view;
- ep.mesh.NbLayer = 1; //this may be more general for defining different layers
- ep.mesh.hLayer.clear();
- ep.mesh.hLayer.push_back(hLayer);
- ep.mesh.NbElmLayer.clear();
- ep.mesh.NbElmLayer.push_back(nbLayers);
- ep.mesh.ExtrudeMesh = true;
- ep.geo.Source = e->tag();
-
+ ExtrudeParams *ep = new ExtrudeParams;
+ ep->mesh.BoundaryLayerIndex = dir;
+ ep->mesh.ViewIndex = view;//view -5 for centerline based extrude
+ ep->mesh.NbLayer = 1; //this may be more general for defining different layers
+ ep->mesh.hLayer.clear();
+ ep->mesh.hLayer.push_back(hLayer);
+ ep->mesh.NbElmLayer.clear();
+ ep->mesh.NbElmLayer.push_back(nbLayers);
+ ep->mesh.ExtrudeMesh = true;
+ ep->geo.Source = e->tag();
+
int type = BOUNDARY_LAYER;
double T0=0., T1=0., T2=0.;
double A0=0., A1=0., A2=0.;
@@ -313,22 +313,22 @@ GEntity* GeoFactory::extrudeBoundaryLayer(GModel *gm, GEntity *e, int nbLayers,
shape.Type = v->Typ;
}
else if (e->dim() == 1){
- ((GEdge*)e)->meshAttributes.extrude = &ep;
+ ((GEdge*)e)->meshAttributes.extrude = ep;
Curve *c = FindCurve(e->tag());
if(!c) printf("curve %d not found \n", e->tag());
shape.Num = c->Num;
shape.Type = c->Typ;
}
else if (e->dim() == 2){
- ((GFace*)e)->meshAttributes.extrude = &ep;
+ ((GFace*)e)->meshAttributes.extrude = ep;
Surface *s = FindSurface(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));
+ // 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);
@@ -338,7 +338,6 @@ GEntity* GeoFactory::extrudeBoundaryLayer(GModel *gm, GEntity *e, int nbLayers,
printf("surface %d NOT found \n", e->tag());
exit(1);
}
-
shape.Num = s->Num;
shape.Type = s->Typ;
}
@@ -351,7 +350,7 @@ GEntity* GeoFactory::extrudeBoundaryLayer(GModel *gm, GEntity *e, int nbLayers,
T0, T1, T2,
A0, A1, A2,
X0, X1, X2, alpha,
- &ep,
+ ep,
list_out);
//create GEntities
diff --git a/Geo/Geo.cpp b/Geo/Geo.cpp
index 6a8dbdd33eba9ea0ea50e1bffc1abd266e0bcb81..8b77da198a7072b75e179851959813de19ad445c 100644
--- a/Geo/Geo.cpp
+++ b/Geo/Geo.cpp
@@ -2763,6 +2763,7 @@ void ExtrudeShapes(int type, List_T *list_in,
case MSH_SURF_TRIC:
case MSH_SURF_PLAN:
case MSH_SURF_DISCRETE:
+ case MSH_SURF_COMPOUND:
{
Volume *pv = 0;
Shape top;
@@ -2806,9 +2807,10 @@ void ExtrudeShapes(int type, List_T *list_in,
static int compareTwoPoints(const void *a, const void *b)
{
+
Vertex *q = *(Vertex **)a;
Vertex *w = *(Vertex **)b;
-
+
if(q->Typ != w->Typ)
return q->Typ - w->Typ;
@@ -2871,7 +2873,8 @@ static int compareTwoSurfaces(const void *a, const void *b)
// checking types is the "right thing" to do (see e.g. compareTwoCurves)
// but it would break backward compatibility (see e.g. tutorial/t2.geo),
// so let's just do it for boundary layer surfaces for now:
- if(s1->Typ == MSH_SURF_BND_LAYER || s2->Typ == MSH_SURF_BND_LAYER){
+ if(s1->Typ == MSH_SURF_BND_LAYER || s2->Typ == MSH_SURF_BND_LAYER ||
+ s1->Typ == MSH_SURF_COMPOUND || s2->Typ == MSH_SURF_COMPOUND ){
if(s1->Typ != s2->Typ) return s1->Typ - s2->Typ;
}
@@ -3020,6 +3023,7 @@ static void ReplaceDuplicatePoints()
Tree_Action(points2delete, Free_Vertex);
Tree_Delete(points2delete);
Tree_Delete(allNonDuplicatedPoints);
+
}
static void ReplaceDuplicateCurves()
diff --git a/Mesh/CenterlineField.cpp b/Mesh/CenterlineField.cpp
index dcba3231a16d397f2a506a1225dcd408b347e504..9d63f325536c3f72737eb466605c5aac42f4183f 100644
--- a/Mesh/CenterlineField.cpp
+++ b/Mesh/CenterlineField.cpp
@@ -840,8 +840,10 @@ void Centerline::extrudeBoundaryLayerWall(){
for (int i= 0; i< NF; i++){
GFace *gf = current->getFaceByTag(NF+i+1);//at this moment compound is not meshed yet exist yet
+
+ int dir = 0;
current->setFactory("Gmsh");
- current->extrudeBoundaryLayer(gf, nbElemLayer, hLayer, 0, -1);
+ current->extrudeBoundaryLayer(gf, nbElemLayer, hLayer, dir, -1);
//view -5 to scale hLayer by radius in BoundaryLayers.cpp
}
diff --git a/Mesh/Generator.cpp b/Mesh/Generator.cpp
index 578fea9eb65543742b75444deb81527b3bf80980..f89e9a4622cc1204938d9d619c6a0159d046d0fd 100644
--- a/Mesh/Generator.cpp
+++ b/Mesh/Generator.cpp
@@ -29,6 +29,7 @@
#include "Generator.h"
#include "meshGFaceLloyd.h"
#include "CenterlineField.h"
+#include "GFaceCompound.h"
#include "Field.h"
#include "Options.h"
#include "simple3D.h"
@@ -473,7 +474,6 @@ static void Mesh2D(GModel *m)
}
}
-
#if defined(HAVE_BFGS)
// lloyd optimization
if (CTX::instance()->mesh.optimizeLloyd){
diff --git a/Mesh/meshGFace.cpp b/Mesh/meshGFace.cpp
index 12fc4e48e94f2ee4251beaa4ec3c7b06c4406b63..8bff17cb4411f5e51782a49dc89431f5315bfdb9 100644
--- a/Mesh/meshGFace.cpp
+++ b/Mesh/meshGFace.cpp
@@ -2118,8 +2118,11 @@ void orientMeshGFace::operator()(GFace *gf)
if(gf->geomType() == GEntity::DiscreteSurface) return;
if(gf->geomType() == GEntity::ProjectionFace) return;
- if(gf->geomType() == GEntity::CompoundSurface) return;
if(gf->geomType() == GEntity::BoundaryLayerSurface) return;
+ if(gf->geomType() == GEntity::CompoundSurface ) {
+ GFaceCompound *gfc = (GFaceCompound*) gf;
+ if (gfc->getCompounds().size() != 1) return;
+ }
if(!gf->getNumMeshElements()) return;
diff --git a/Mesh/meshGFaceOptimize.cpp b/Mesh/meshGFaceOptimize.cpp
index e5106e97b973d8cbcf38db26228f722d966e60d2..24ef02767e44cd1dc7f85db02373ef8228e43044 100644
--- a/Mesh/meshGFaceOptimize.cpp
+++ b/Mesh/meshGFaceOptimize.cpp
@@ -1054,7 +1054,7 @@ struct quadBlob {
v2t_cont :: const_iterator it = adj.find(v);
int ss = temp.size();
std::vector<MElement*> elems = it->second;
- /*
+
//EMI: try to orient BNodes the same as quad orientations
for (int j=0;j<elems.size();j++){
bool found =false;
@@ -1072,33 +1072,8 @@ struct quadBlob {
}
if (found) break;
}
- */
- for (int j=0;j<elems.size();j++){
- bool found = false;
- for (int i=0;i<4;i++){
- MEdge e = it->second[j]->getEdge(i);
- if (e.getVertex(0) == v &&
- inBoundary(e.getVertex(1),bnodes) &&
- !inBoundary(e.getVertex(1),temp)) {
- v = e.getVertex(1);
- temp.push_back(e.getVertex(1));
- found = true;
- break;
- }
- else if (e.getVertex(1) == v &&
- inBoundary(e.getVertex(0),bnodes) &&
- !inBoundary(e.getVertex(0),temp)) {
- v = e.getVertex(0);
- temp.push_back(e.getVertex(0));
- found = true;
- break;
- }
- }
- if (found)break;
- }
-
- //JF this does not orient
+ //JF this does not orient quads
// for (int j=0;j<elems.size();j++){
// bool found = false;
// if(inBlob(elems[j])){
@@ -1112,14 +1087,14 @@ struct quadBlob {
// found = true;
// break;
// }
- // //else if (e.getVertex(1) == v &&
- // // inBoundary(e.getVertex(0),bnodes) &&
- // // !inBoundary(e.getVertex(0),temp)) {
- // // v = e.getVertex(0);
- // // temp.push_back(e.getVertex(0));
- // // found = true;
- // // break;
- // // }
+ // else if (e.getVertex(1) == v &&
+ // inBoundary(e.getVertex(0),bnodes) &&
+ // !inBoundary(e.getVertex(0),temp)) {
+ // v = e.getVertex(0);
+ // temp.push_back(e.getVertex(0));
+ // found = true;
+ // break;
+ // }
// }
// }
// if (found)break;
diff --git a/benchmarks/3d/percolation.geo b/benchmarks/3d/percolation.geo
index 836624cde804cedfb6d3d88165e6eaa65c352493..e2c197a5f631c499483441d0b444b622d708012c 100644
--- a/benchmarks/3d/percolation.geo
+++ b/benchmarks/3d/percolation.geo
@@ -133,3 +133,4 @@ c = newsl; Surface Loop(c) = {43,-1,35,39,48,47};
v = newv; Volume(v) = {c, holes[]};
Coherence; // make sure duplicate nodes are removed when sphere_dist=0
+Field[1] = BoundaryLayer;