diff --git a/Geo/GFaceCompound.cpp b/Geo/GFaceCompound.cpp
index c2a9e97c84a8d4e536873a09c2a2aeb7fd2ad978..dbba6933626166a994b498e1efb6b33b07f09297 100644
--- a/Geo/GFaceCompound.cpp
+++ b/Geo/GFaceCompound.cpp
@@ -995,7 +995,6 @@ bool GFaceCompound::parametrize() const
}
printStuff(55);
oriented = false; //checkOrientation(0);
- printStuff(66);
if (!oriented) oriented = checkOrientation(0, true);
printStuff(77);
if (_type==SPECTRAL && (!oriented || checkOverlap(vert)) ){
@@ -1644,11 +1643,11 @@ bool GFaceCompound::parametrize_conformal_spectral() const
myAssembler.numberVertex(v, 0, 2);
}
- for(std::set<MVertex *>::iterator itv = fillNodes.begin(); itv !=fillNodes.end() ; ++itv){
- MVertex *v = *itv;
- myAssembler.numberVertex(v, 0, 1);
- myAssembler.numberVertex(v, 0, 2);
- }
+ // for(std::set<MVertex *>::iterator itv = fillNodes.begin(); itv !=fillNodes.end() ; ++itv){
+ // MVertex *v = *itv;
+ // myAssembler.numberVertex(v, 0, 1);
+ // myAssembler.numberVertex(v, 0, 2);
+ // }
laplaceTerm laplace1(model(), 1, ONE);
laplaceTerm laplace2(model(), 2, ONE);
@@ -1665,13 +1664,13 @@ bool GFaceCompound::parametrize_conformal_spectral() const
}
}
- for (std::list<MTriangle*>::iterator it2 = fillTris.begin(); it2 != fillTris.end(); it2++){
- SElement se((*it2));
- laplace1.addToMatrix(myAssembler, &se);
- laplace2.addToMatrix(myAssembler, &se);
- cross12.addToMatrix(myAssembler, &se);
- cross21.addToMatrix(myAssembler, &se);
- }
+ // for (std::list<MTriangle*>::iterator it2 = fillTris.begin(); it2 != fillTris.end(); it2++){
+ // SElement se((*it2));
+ // laplace1.addToMatrix(myAssembler, &se);
+ // laplace2.addToMatrix(myAssembler, &se);
+ // cross12.addToMatrix(myAssembler, &se);
+ // cross21.addToMatrix(myAssembler, &se);
+ // }
double epsilon = 1.e-6;
for(std::set<MVertex *>::iterator itv = allNodes.begin(); itv !=allNodes.end() ; ++itv){
@@ -2581,7 +2580,7 @@ double GFaceCompound::checkAspectRatio() const
}
void GFaceCompound::coherencePatches() const
-{
+{
if (_mapping == RBF) return;
Msg::Info("Re-orient all %d compound patches normals coherently",
_compound.size());
@@ -2612,31 +2611,31 @@ void GFaceCompound::coherencePatches() const
}
std::set<MElement* , std::less<MElement*> > touched;
-int iE, si, iE2, si2;
-touched.insert(allElems[0]);
-while(touched.size() != allElems.size()){
- for(unsigned int i = 0; i < allElems.size(); i++){
- MElement *t = allElems[i];
- std::set<MElement*, std::less<MElement*> >::iterator it2 = touched.find(t);
- if(it2 != touched.end()){
- for (int j = 0; j < t->getNumEdges(); j++){
- MEdge me = t->getEdge(j);
- t->getEdgeInfo(me, iE,si);
- std::map<MEdge, std::set<MElement*>, Less_Edge >::iterator it =
- edge2elems.find(me);
- std::set<MElement*, std::less<MElement*> > mySet = it->second;
- for(std::set<MElement*, std::less<MElement*> >::iterator itt = mySet.begin();
- itt != mySet.end(); itt++){
- if (*itt != t){
- (*itt)->getEdgeInfo(me,iE2,si2);
- if(si == si2) (*itt)->revert();
- touched.insert(*itt);
+ int iE, si, iE2, si2;
+ touched.insert(allElems[0]);
+ while(touched.size() != allElems.size()){
+ for(unsigned int i = 0; i < allElems.size(); i++){
+ MElement *t = allElems[i];
+ std::set<MElement*, std::less<MElement*> >::iterator it2 = touched.find(t);
+ if(it2 != touched.end()){
+ for (int j = 0; j < t->getNumEdges(); j++){
+ MEdge me = t->getEdge(j);
+ t->getEdgeInfo(me, iE,si);
+ std::map<MEdge, std::set<MElement*>, Less_Edge >::iterator it =
+ edge2elems.find(me);
+ std::set<MElement*, std::less<MElement*> > mySet = it->second;
+ for(std::set<MElement*, std::less<MElement*> >::iterator itt = mySet.begin();
+ itt != mySet.end(); itt++){
+ if (*itt != t){
+ (*itt)->getEdgeInfo(me,iE2,si2);
+ if(si == si2) { (*itt)->revert();}
+ touched.insert(*itt);
+ }
}
}
}
}
}
- }
}
void GFaceCompound::coherenceNormals()
diff --git a/Geo/GFaceCompound.h b/Geo/GFaceCompound.h
index dad482fd9c72de3075a2a9e3f766bd82a98cc1cc..93ae99228e3b0ad1af9f0518247fa29b58fafd31 100644
--- a/Geo/GFaceCompound.h
+++ b/Geo/GFaceCompound.h
@@ -126,8 +126,7 @@ class GFaceCompound : public GFace {
void getBoundingEdges();
void getUniqueEdges(std::set<GEdge*> &_unique);
void computeALoop(std::set<GEdge*> &_unique, std::list<GEdge*> &);
- void getTriangle(double u, double v, GFaceCompoundTriangle **lt,
- double &_u, double &_v) const;
+
virtual double locCurvature(MTriangle *t, double u, double v) const;
double getSizeH() const;
@@ -160,6 +159,8 @@ class GFaceCompound : public GFace {
virtual GEntity::GeomType geomType() const { return CompoundSurface; }
ModelType getNativeType() const { return GmshModel; }
void * getNativePtr() const { return 0; }
+ void getTriangle(double u, double v, GFaceCompoundTriangle **lt,
+ double &_u, double &_v) const;
virtual SPoint2 getCoordinates(MVertex *v) const;
virtual double curvatureMax(const SPoint2 ¶m) const;
virtual double curvatures(const SPoint2 ¶m, SVector3 *dirMax, SVector3 *dirMin,
diff --git a/Mesh/CenterlineField.cpp b/Mesh/CenterlineField.cpp
index 9d63f325536c3f72737eb466605c5aac42f4183f..90b36f2f52a557b99be408001c43c33f07324aa6 100644
--- a/Mesh/CenterlineField.cpp
+++ b/Mesh/CenterlineField.cpp
@@ -647,7 +647,7 @@ void Centerline::createSplitCompounds(){
Msg::Info("Create Compound Surface (%d) = %d discrete face",
num_gfc, pf->tag());
- GFace *gfc = current->addCompoundFace(f_compound, 1, 0, num_gfc); //1=conf_spectral 4=convex_circle
+ GFace *gfc = current->addCompoundFace(f_compound, 7, 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;
@@ -841,7 +841,14 @@ 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;
+ 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);
//view -5 to scale hLayer by radius in BoundaryLayers.cpp
diff --git a/Mesh/meshGFace.cpp b/Mesh/meshGFace.cpp
index 8bff17cb4411f5e51782a49dc89431f5315bfdb9..0e328df6492dffdf29a2e3da05dcd1b3dbc72bd2 100644
--- a/Mesh/meshGFace.cpp
+++ b/Mesh/meshGFace.cpp
@@ -17,6 +17,7 @@
#include "GVertex.h"
#include "GEdge.h"
#include "GEdgeCompound.h"
+#include "robustPredicates.h"
#include "GFace.h"
#include "GModel.h"
#include "MVertex.h"
@@ -1787,7 +1788,7 @@ void deMeshGFace::operator() (GFace *gf)
}
// for debugging, change value from -1 to -100;
-int debugSurface = -1; //-1;
+int debugSurface = -100; //-1;
void meshGFace::operator() (GFace *gf)
{
@@ -2119,13 +2120,46 @@ void orientMeshGFace::operator()(GFace *gf)
if(gf->geomType() == GEntity::DiscreteSurface) return;
if(gf->geomType() == GEntity::ProjectionFace) 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;
+ 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());}
+
+ // 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;
// we could add the following to provide backward compatibility:
// if(gf->meshAttributes.Method == MESH_TRANSFINITE) return;
@@ -2181,7 +2215,7 @@ void orientMeshGFace::operator()(GFace *gf)
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());
+ Msg::Debug("Reverting 2 orientation of mesh in face %d", gf->tag());
for(unsigned int k = 0; k < gf->getNumMeshElements(); k++)
gf->getMeshElement(k)->revert();
}
diff --git a/Solver/multiscaleLaplace.cpp b/Solver/multiscaleLaplace.cpp
index a0a5c9b2a6d9cb62850009416f79a8e3c3f981ab..aaaf4ec24ee406f73d5a85ac82acd6a6519f54b1 100644
--- a/Solver/multiscaleLaplace.cpp
+++ b/Solver/multiscaleLaplace.cpp
@@ -876,8 +876,8 @@ multiscaleLaplace::multiscaleLaplace (std::vector<MElement *> &elements,
//Split the mesh in left and right
//or Cut the mesh in left and right
- //splitElems(elements);
- cutElems(elements);
+ splitElems(elements);
+ //cutElems(elements);
}