diff --git a/Geo/GModel.cpp b/Geo/GModel.cpp
index 5c522ef45c6c39d473bbc99799f9bc2851b94a2a..658660c04a13ebe5347228c8e29081517dc54726 100644
--- a/Geo/GModel.cpp
+++ b/Geo/GModel.cpp
@@ -580,7 +580,7 @@ int GModel::adaptMesh(std::vector<int> technique, std::vector<simpleFunction<dou
//adapt only upper most dimension
else{
- while(1) {
+ while(1) {
Msg::Info("-- adaptMesh ITER =%d ", ITER);
std::vector<MElement*> elements;
diff --git a/Mesh/BoundaryLayers.cpp b/Mesh/BoundaryLayers.cpp
index 02c6b467b555f94c4b0622a89221f44a85cc8468..85a9cd350be968dd386dbc9e0ba76a73a7ac0f40 100644
--- a/Mesh/BoundaryLayers.cpp
+++ b/Mesh/BoundaryLayers.cpp
@@ -12,6 +12,7 @@
#include "meshGEdge.h"
#include "meshGFace.h"
#include "GmshMessage.h"
+#include "Field.h"
#if defined(HAVE_POST)
#include "PView.h"
@@ -68,7 +69,7 @@ template<class T>
static void addExtrudeNormals(std::set<T*> &entities,
std::map<int, infoset> &infos)
{
- bool normalize = true, special3dbox = false;
+ bool normalize = true, special3dbox = false, extrudeField=false;
std::vector<OctreePost*> octrees;
for(typename std::set<T*>::iterator it = entities.begin(); it != entities.end(); it++){
@@ -94,6 +95,10 @@ static void addExtrudeNormals(std::set<T*> &entities,
// build nice 3D "boxes")
special3dbox = true;
}
+ else if(view == -5){
+ // Force extrusion normals with a field
+ extrudeField = true;
+ }
else
Msg::Error("Unknown View[%d]: using normals instead", view);
}
@@ -133,6 +138,18 @@ static void addExtrudeNormals(std::set<T*> &entities,
}
}
}
+ if(extrudeField){ // force normals according to field
+ for(smooth_data::iter it = ExtrudeParams::normals[i]->begin();
+ it != ExtrudeParams::normals[i]->end(); it++){
+ GEntity *ge = (GEntity*)(*entities.begin());
+ FieldManager *fields = ge->model()->getFields();
+ if(fields->getBackgroundField() > 0){
+ Field *f = fields->get(fields->getBackgroundField());
+ double radius = (*f)(it->x,it->y,it->z);
+ for(int k = 0; k < 3; k++) it->vals[k] *= radius;
+ }
+ }
+ }
#if defined(HAVE_POST)
if(octrees.size()){ // scale normals by scalar views
for(smooth_data::iter it = ExtrudeParams::normals[i]->begin();
diff --git a/Mesh/CenterlineField.cpp b/Mesh/CenterlineField.cpp
index 90b36f2f52a557b99be408001c43c33f07324aa6..5b14a108877b6213e0b77bc9e111b2e5ef77941d 100644
--- a/Mesh/CenterlineField.cpp
+++ b/Mesh/CenterlineField.cpp
@@ -838,19 +838,22 @@ void Centerline::extrudeBoundaryLayerWall(){
Msg::Info("Centerline: extrude boundary layer wall %d %g ", nbElemLayer, hLayer);
+ //orient extrude direction outward
+ int dir = 0;
+ MElement *e = current->getFaceByTag(1)->getMeshElement(0);
+ SVector3 ne = e->getFace(0).normal();
+ SVector3 ps(e->getVertex(0)->x(), e->getVertex(0)->y(), e->getVertex(0)->z());
+ double xyz[3] = {ps.x(), ps.y(), ps.z()};
+ kdtree->annkSearch(xyz, 1, index, dist);
+ SVector3 pc(nodes[index[0]][0], nodes[index[0]][1], nodes[index[0]][2]);
+ SVector3 nc = ps-pc;
+ if (dot(ne,nc) < 0) dir = 1;
+ if (dir ==1 && hLayer > 0 ) hLayer *= -1.0;
+
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 = 1;
-
- //do sthg here
- //MElement *e = current->getFaceByTag(i+1)->getMeshElement(i);
-
-
- if (dir ==1 && hLayer > 0 ) hLayer *= -1.0;
- printf("dir=%d hlayer =%g \n", dir, hLayer);
current->setFactory("Gmsh");
- current->extrudeBoundaryLayer(gf, nbElemLayer, hLayer, dir, -1);
+ current->extrudeBoundaryLayer(gf, nbElemLayer, hLayer, dir, -5);
//view -5 to scale hLayer by radius in BoundaryLayers.cpp
}
@@ -1078,26 +1081,32 @@ double Centerline::operator() (double x, double y, double z, GEntity *ge){
double xyz[3] = {x,y,z};
//take xyz = closest point on boundary in case we are on the planar in/out faces
- bool isCompound = (ge->dim() == 2 && ge->geomType() == GEntity::CompoundSurface) ? true: false;
- std::list<GFace*> cFaces;
- if (isCompound) cFaces = ((GFaceCompound*)ge)->getCompounds();
- if ( ge->dim() == 3 || (ge->dim() == 2 && ge->geomType() == GEntity::Plane) ||
- (isCompound && (*cFaces.begin())->geomType() == GEntity::Plane) ){
- int num_neighbours = 1;
- kdtreeR->annkSearch(xyz, num_neighbours, index, dist);
- xyz[0] = nodesR[index[0]][0];
- xyz[1] = nodesR[index[0]][1];
- xyz[2] = nodesR[index[0]][2];
+ bool isCompound = false;
+ if(ge){
+ if (ge->dim() == 2 && ge->geomType() == GEntity::CompoundSurface) isCompound = true;
+ std::list<GFace*> cFaces;
+ if (isCompound) cFaces = ((GFaceCompound*)ge)->getCompounds();
+ if ( ge->dim() == 3 || (ge->dim() == 2 && ge->geomType() == GEntity::Plane) ||
+ (isCompound && (*cFaces.begin())->geomType() == GEntity::Plane) ){
+ int num_neighbours = 1;
+ kdtreeR->annkSearch(xyz, num_neighbours, index, dist);
+ xyz[0] = nodesR[index[0]][0];
+ xyz[1] = nodesR[index[0]][1];
+ xyz[2] = nodesR[index[0]][2];
+ }
}
int num_neighbours = 1;
kdtree->annkSearch(xyz, num_neighbours, index, dist);
double d = sqrt(dist[0]);
double lc = 2*M_PI*d/nbPoints;
- return lc;
+
+ if(!ge) { return d;}
+ else return lc;
}
+
void Centerline::operator() (double x, double y, double z, SMetric3 &metr, GEntity *ge){
if (update_needed){
diff --git a/Mesh/CenterlineField.h b/Mesh/CenterlineField.h
index a8adc32a724742bd291b14cf48cd4411f115383a..e39806bbaacffd792c9c1236731b8326c2792d48 100644
--- a/Mesh/CenterlineField.h
+++ b/Mesh/CenterlineField.h
@@ -117,7 +117,7 @@ class Centerline : public Field{
double operator() (double x, double y, double z, GEntity *ge=0);
//anisotropic operator
void operator() (double x, double y, double z, SMetric3 &metr, GEntity *ge=0);
-
+
//temporary operator where v1, v2 and v3 are three orthonormal directions
void operator()(double x,double y,double z,SVector3& v1,SVector3& v2,SVector3& v3,GEntity* ge=0);
diff --git a/Mesh/meshGFace.cpp b/Mesh/meshGFace.cpp
index 6d88dd5090e5e4498791ba89748b8b5feb9ae85a..99bf99e2ab87a7de24a0817bd685338a11df3cd7 100644
--- a/Mesh/meshGFace.cpp
+++ b/Mesh/meshGFace.cpp
@@ -2123,41 +2123,41 @@ void orientMeshGFace::operator()(GFace *gf)
if(!gf->getNumMeshElements()) return;
- if(gf->geomType() == GEntity::CompoundSurface ) return;
+ //if(gf->geomType() == GEntity::CompoundSurface ) return;
//do sthg for compound face
- // if(gf->geomType() == GEntity::CompoundSurface ) {
- // GFaceCompound *gfc = (GFaceCompound*) gf;
- // //if (gfc->getCompounds().size() != 1) return;
- // //else{printf("compound face %d orient \n", gfc->tag());}
+ if(gf->geomType() == GEntity::CompoundSurface ) {
+ GFaceCompound *gfc = (GFaceCompound*) gf;
+ //if (gfc->getCompounds().size() != 1) return;
+ //else{printf("compound face %d orient \n", gfc->tag());}
- // std::list<GFace*> comp = gfc->getCompounds();
- // MTriangle *lt = (*comp.begin())->triangles[0];
- // SPoint2 c0 = gfc->getCoordinates(lt->getVertex(0));
- // SPoint2 c1 = gfc->getCoordinates(lt->getVertex(1));
- // SPoint2 c2 = gfc->getCoordinates(lt->getVertex(2));
- // double p0[2] = {c0[0],c0[1]};
- // double p1[2] = {c1[0],c1[1]};
- // double p2[2] = {c2[0],c2[1]};
- // double normal = robustPredicates::orient2d(p0, p1, p2);
-
- // MElement *e = gfc->getMeshElement(0);
- // SPoint2 v1,v2,v3;
- // reparamMeshVertexOnFace(e->getVertex(0), gf, v1, false);
- // reparamMeshVertexOnFace(e->getVertex(1), gf, v2, false);
- // reparamMeshVertexOnFace(e->getVertex(2), gf, v3, false);
- // SVector3 C1(v1.x(), v1.y(), 0.0);
- // SVector3 C2(v2.x(), v2.y(), 0.0);
- // SVector3 C3(v3.x(), v3.y(), 0.0);
- // SVector3 n1 = crossprod(C2-C1,C3-C1);
-
- // if(normal*n1.z() < 0){
- // Msg::Debug("Reverting orientation of mesh in compound face %d", gf->tag());
- // for(unsigned int k = 0; k < gf->getNumMeshElements(); k++)
- // gfc->getMeshElement(k)->revert();
- // }
- // return;
-
- // }
+ std::list<GFace*> comp = gfc->getCompounds();
+ MTriangle *lt = (*comp.begin())->triangles[0];
+ SPoint2 c0 = gfc->getCoordinates(lt->getVertex(0));
+ SPoint2 c1 = gfc->getCoordinates(lt->getVertex(1));
+ SPoint2 c2 = gfc->getCoordinates(lt->getVertex(2));
+ double p0[2] = {c0[0],c0[1]};
+ double p1[2] = {c1[0],c1[1]};
+ double p2[2] = {c2[0],c2[1]};
+ double normal = robustPredicates::orient2d(p0, p1, p2);
+
+ MElement *e = gfc->getMeshElement(0);
+ SPoint2 v1,v2,v3;
+ reparamMeshVertexOnFace(e->getVertex(0), gf, v1, false);
+ reparamMeshVertexOnFace(e->getVertex(1), gf, v2, false);
+ reparamMeshVertexOnFace(e->getVertex(2), gf, v3, false);
+ SVector3 C1(v1.x(), v1.y(), 0.0);
+ SVector3 C2(v2.x(), v2.y(), 0.0);
+ SVector3 C3(v3.x(), v3.y(), 0.0);
+ SVector3 n1 = crossprod(C2-C1,C3-C1);
+
+ if(normal*n1.z() < 0){
+ Msg::Debug("Reverting orientation of mesh in compound face %d", gf->tag());
+ for(unsigned int k = 0; k < gf->getNumMeshElements(); k++)
+ gfc->getMeshElement(k)->revert();
+ }
+ return;
+
+ }
// In old versions we did not reorient transfinite surface meshes;
diff --git a/Mesh/meshMetric.cpp b/Mesh/meshMetric.cpp
index 8ba87063ad9318234fa6bc0e152785ca955b3d40..242d3a52e7a94c7e63b0b3af4778b4e2dd2c3124 100644
--- a/Mesh/meshMetric.cpp
+++ b/Mesh/meshMetric.cpp
@@ -371,8 +371,8 @@ void meshMetric::computeMetric(){
//printf("%d elements are considered in the meshMetric \n",(int)_elements.size());
computeValues();
- computeHessian_FE();
- //computeHessian_LS();
+ //computeHessian_FE();
+ computeHessian_LS();
int metricNumber = setOfMetrics.size();