diff --git a/Common/OS.cpp b/Common/OS.cpp
index fda2a9afcdb1428e49d4dbbc9b3ca1868e76643e..4b21dff83990927ddb9bde8e88575f4ed2bf1cda 100644
--- a/Common/OS.cpp
+++ b/Common/OS.cpp
@@ -348,12 +348,12 @@ void CheckResources()
// Try to get at least 16 MB of stack. Running with too small a stack
// can cause crashes in the recursive calls (e.g. for tet
// classification in 3D Delaunay)
- if(r.rlim_cur < 16 * 1024 * 1024) {
- Msg::Info("Increasing process stack size (%d kB < 16 MB)",
- r.rlim_cur / 1024);
- r.rlim_cur = r.rlim_max;
- setrlimit(RLIMIT_STACK, &r);
- }
+ if(r.rlim_cur < 16 * 1024 * 1024) {
+ Msg::Info("Increasing process stack size (%d kB < 16 MB)",
+ r.rlim_cur / 1024);
+ r.rlim_cur = r.rlim_max;
+ setrlimit(RLIMIT_STACK, &r);
+ }
#endif
}
diff --git a/Fltk/classificationEditor.cpp b/Fltk/classificationEditor.cpp
index 63931d77b45a9a27aa4dd3466e417d46cd6c2074..6088536cc51537b1cc06f034a1f38be35c980c75 100644
--- a/Fltk/classificationEditor.cpp
+++ b/Fltk/classificationEditor.cpp
@@ -262,6 +262,7 @@ static void reset_selection_cb(Fl_Widget *w, void *data)
static void classify_cb(Fl_Widget *w, void *data)
{
+ // printf("cpoc<<<\n");
classificationEditor *e = (classificationEditor *)data;
if(!e->selected) {
@@ -270,9 +271,14 @@ static void classify_cb(Fl_Widget *w, void *data)
GModel::current()->getMaxElementaryNumber(1) + 1, 0, 0);
GModel::current()->add(e->selected);
}
+ // printf("cpoc<<<\n");
- if(e->toggles[CLASS_TOGGLE_ENSURE_PARAMETRIZABLE_SURFACES]->value())
+ std::map<MVertex* , std::pair<SVector3, SVector3> > CURVATURES;
+ if(e->toggles[CLASS_TOGGLE_ENSURE_PARAMETRIZABLE_SURFACES]->value()){
+ computeDiscreteCurvatures(GModel::current(), CURVATURES );
computeEdgeCut(GModel::current(), e->selected->lines, 100000);
+ }
+ // printf("cpoc<<<\n");
computeNonManifoldEdges(GModel::current(), e->selected->lines, true);
@@ -285,16 +291,21 @@ static void classify_cb(Fl_Widget *w, void *data)
delete e->selected;
e->selected = 0;
}
+ // printf("cpoc<<<\n");
e->elements.clear();
e->edges_detected.clear();
GModel::current()->pruneMeshVertexAssociations();
NoElementsSelectedMode(e);
+ // printf("cpoc<<< PARAMETRIZE\n");
if(e->toggles[CLASS_TOGGLE_ENSURE_PARAMETRIZABLE_SURFACES]->value()){
+ // printf("ZAZOU\n");
parametrizeAllGEdge(GModel::current());
- parametrizeAllGFace(GModel::current());
+ // printf("ZAZOU\n");
+ parametrizeAllGFace(GModel::current(), &CURVATURES);
}
+ // printf("cpoc<<<\n");
}
diff --git a/Geo/GFace.h b/Geo/GFace.h
index 11c14d2ed0963fd1b673c2b8d9ce86f495a8c186..f1ca7c23220da3a6ceb2d1266d82b1f7a514947e 100644
--- a/Geo/GFace.h
+++ b/Geo/GFace.h
@@ -327,6 +327,7 @@ public:
std::vector<SPoint2> stl_vertices_uv;
std::vector<SPoint3> stl_vertices_xyz;
std::vector<SVector3> stl_normals;
+ std::vector<SVector3> stl_curvatures;
std::vector<int> stl_triangles;
// a vertex array containing a geometrical representation of the
diff --git a/Geo/GModel.cpp b/Geo/GModel.cpp
index e1eff1dc1a4aab8e082991cb828cea470f383964..eba8bb620e78bc59488d80d83be99b0bc8fa6aad 100644
--- a/Geo/GModel.cpp
+++ b/Geo/GModel.cpp
@@ -2938,53 +2938,71 @@ void GModel::classifyAllFaces(double angleThreshold, bool includeBoundary)
}
#if defined(HAVE_MESH)
-void recurClassifyEdges(MTri3 *t, std::map<MTriangle *, GFace *> &reverse,
+static void recurClassifyEdges(MTri3 *t, std::map<MTriangle *, GFace *> &reverse,
std::map<MLine *, GEdge *, compareMLinePtr> &lines,
std::set<MLine *> &touched,
std::set<MTri3 *> &trisTouched,
std::map<std::pair<int, int>, GEdge *> &newEdges)
{
- if(!t->isDeleted()) {
- trisTouched.erase(t);
- t->setDeleted(true);
- GFace *gf1 = reverse[t->tri()];
- for(int i = 0; i < 3; i++) {
- GFace *gf2 = 0;
- MTri3 *tn = t->getNeigh(i);
- if(tn) gf2 = reverse[tn->tri()];
- edgeXface exf(t, i);
- MLine ml(exf._v(0), exf._v(1));
- std::map<MLine *, GEdge *, compareMLinePtr>::iterator it =
- lines.find(&ml);
- if(it != lines.end()) {
- if(touched.find(it->first) == touched.end()) {
- GEdge *ge = getNewModelEdge(gf1, gf2, newEdges);
- if(ge) ge->lines.push_back(it->first);
- touched.insert(it->first);
- }
- }
- if(tn)
- recurClassifyEdges(tn, reverse, lines, touched, trisTouched, newEdges);
- }
- }
-}
-void recurClassify(MTri3 *t, GFace *gf,
+ std::stack<MTri3 *> _stack;
+ _stack.push(t);
+
+ while (!_stack.empty()){
+ t = _stack.top();
+ _stack.pop();
+ // printf("%d\n",rec);
+ if(!t->isDeleted()) {
+ trisTouched.erase(t);
+ t->setDeleted(true);
+ GFace *gf1 = reverse[t->tri()];
+ for(int i = 0; i < 3; i++) {
+ GFace *gf2 = 0;
+ MTri3 *tn = t->getNeigh(i);
+ if(tn) gf2 = reverse[tn->tri()];
+ edgeXface exf(t, i);
+ MLine ml(exf._v(0), exf._v(1));
+ std::map<MLine *, GEdge *, compareMLinePtr>::iterator it =
+ lines.find(&ml);
+ if(it != lines.end()) {
+ if(touched.find(it->first) == touched.end()) {
+ GEdge *ge = getNewModelEdge(gf1, gf2, newEdges);
+ if(ge) ge->lines.push_back(it->first);
+ touched.insert(it->first);
+ }
+ }
+ if(tn){
+ _stack.push(tn);
+ // recurClassifyEdges(tn, reverse, lines, touched, trisTouched, newEdges,rec+1);
+ }
+ }
+ }
+ }
+}
+
+static void recurClassify(MTri3 *t, GFace *gf,
std::map<MLine *, GEdge *, compareMLinePtr> &lines,
std::map<MTriangle *, GFace *> &reverse)
{
- if(!t->isDeleted()) {
- gf->triangles.push_back(t->tri());
- reverse[t->tri()] = gf;
- t->setDeleted(true);
- for(int i = 0; i < 3; i++) {
- MTri3 *tn = t->getNeigh(i);
- if(tn) {
- edgeXface exf(t, i);
- MLine ml(exf._v(0), exf._v(1));
- std::map<MLine *, GEdge *, compareMLinePtr>::iterator it =
- lines.find(&ml);
- if(it == lines.end()) recurClassify(tn, gf, lines, reverse);
+ std::stack<MTri3 *> _stack;
+ _stack.push(t);
+
+ while (!_stack.empty()){
+ t = _stack.top();
+ _stack.pop();
+ if(!t->isDeleted()) {
+ gf->triangles.push_back(t->tri());
+ reverse[t->tri()] = gf;
+ t->setDeleted(true);
+ for(int i = 0; i < 3; i++) {
+ MTri3 *tn = t->getNeigh(i);
+ if(tn) {
+ edgeXface exf(t, i);
+ MLine ml(exf._v(0), exf._v(1));
+ std::map<MLine *, GEdge *, compareMLinePtr>::iterator it =
+ lines.find(&ml);
+ if(it == lines.end()) _stack.push(tn);//recurClassify(tn, gf, lines, reverse);
+ }
}
}
}
@@ -3013,6 +3031,7 @@ void GModel::classifyFaces()
for(size_t i = 0; i < edgesToRemove.size(); ++i) {
remove(edgesToRemove[i]);
}
+ // printf("cpoc<<<s\n");
// create triangle 2 triangle connections
std::map<MTriangle *, GFace *> reverse_old;
@@ -3031,7 +3050,9 @@ void GModel::classifyFaces()
}
}
if(tris.empty()) return;
+ // printf("cpoc<<<ss\n");
connectTriangles(tris);
+ // printf("cpoc<<<sq\n");
// classify
std::map<MTriangle *, GFace *> reverse;
@@ -3053,6 +3074,7 @@ void GModel::classifyFaces()
}
++it;
}
+ // printf("cpoc<<<s\n");
// now we have all faces coloured. If some regions were existing, replace
// their faces by the new ones
@@ -3071,6 +3093,7 @@ void GModel::classifyFaces()
}
(*rit)->set(std::vector<GFace *>(_newFaces.begin(), _newFaces.end()));
}
+ // printf("cpoc<<<s\n");
// color some lines
it = tris.begin();
@@ -3090,6 +3113,7 @@ void GModel::classifyFaces()
recurClassifyEdges(*trisTouched.begin(), reverse, lines, touched,
trisTouched, newEdges);
+ // printf("cpoc<<<s\n");
std::map<discreteFace *, std::vector<int> > newFaceTopology;
// check if new edges should not be splitted
@@ -3186,12 +3210,14 @@ void GModel::classifyFaces()
if(gf2) newFaceTopology[gf2].push_back(newGe->tag());
}
}
+ // printf("edededededecpocss<<<s\n");
std::map<discreteFace *, std::vector<int> >::iterator itFT =
newFaceTopology.begin();
for(; itFT != newFaceTopology.end(); ++itFT) {
itFT->first->setBoundEdges(itFT->second);
}
+ // printf("edededededecpocss<<<s\n");
for(std::map<std::pair<int, int>, GEdge *>::iterator it = newEdges.begin();
it != newEdges.end(); ++it) {
@@ -3206,6 +3232,7 @@ void GModel::classifyFaces()
++it;
}
+ // printf("cpoc<<<s\n");
// delete empty mesh faces and reclasssify
std::set<GFace *, GEntityLessThan> fac = faces;
for(fiter fit = fac.begin(); fit != fac.end(); ++fit) {
diff --git a/Geo/GModelIO_MSH4.cpp b/Geo/GModelIO_MSH4.cpp
index fb08c8b9b617b66d4f6b5458c8f33ec242895039..599731f35123412b85a074067cad2112ab140c6f 100644
--- a/Geo/GModelIO_MSH4.cpp
+++ b/Geo/GModelIO_MSH4.cpp
@@ -1297,11 +1297,13 @@ static bool readMSH4Parametrizations(GModel *const model, FILE *fp, bool binary
gf->stl_vertices_uv.clear();
gf->stl_normals.clear();
for (int i=0;i<n;i++){
- double u,v,x,y,z;
- fscanf (fp,"%lf %lf %lf %lf %lf",&u,&v,&x,&y,&z);
+ double u,v,x,y,z,cxM,cyM,czM,cxm,cym,czm;
+ fscanf (fp,"%lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf",&u,&v,&x,&y,&z,&cxM,&cyM,&czM,&cxm,&cym,&czm);
gf->stl_vertices_uv.push_back(SPoint2(u,v));
gf->stl_vertices_xyz.push_back(SPoint3(x,y,z));
gf->stl_normals.push_back(SVector3(0,0,0));
+ gf->stl_curvatures.push_back(SVector3(cxM,cyM,czM));
+ gf->stl_curvatures.push_back(SVector3(cxm,cym,czm));
}
gf->stl_triangles.clear();
for (int i=0;i<t;i++){
@@ -2821,8 +2823,15 @@ static void writeMSH4Parametrizations(GModel *const model, FILE *fp, bool binar
if (gf->stl_vertices_uv.size()){
fprintf(fp,"%d %lu %lu\n",gf->tag(),gf->stl_vertices_uv.size(),gf->stl_triangles.size()/3);
for (size_t i=0;i<gf->stl_vertices_uv.size();i++){
- fprintf(fp,"%g %g %g %g %g\n",gf->stl_vertices_uv[i].x(),gf->stl_vertices_uv[i].y(),
- gf->stl_vertices_xyz[i].x(),gf->stl_vertices_xyz[i].y(),gf->stl_vertices_xyz[i].z());
+ if (gf->stl_curvatures.empty())
+ fprintf(fp,"%22.15e %22.15e %22.15e %22.15e %22.15e 0 0 0 0 0 0\n",gf->stl_vertices_uv[i].x(),gf->stl_vertices_uv[i].y(),
+ gf->stl_vertices_xyz[i].x(),gf->stl_vertices_xyz[i].y(),gf->stl_vertices_xyz[i].z());
+ else
+ fprintf(fp,"%22.15e %22.15e %22.15e %22.15e %22.15e %22.15e %22.15e %22.15e %22.15e %22.15e %22.15e\n",gf->stl_vertices_uv[i].x(),gf->stl_vertices_uv[i].y(),
+ gf->stl_vertices_xyz[i].x(),gf->stl_vertices_xyz[i].y(),gf->stl_vertices_xyz[i].z(),
+ gf->stl_curvatures[2*i].x(),gf->stl_curvatures[2*i].y(),gf->stl_curvatures[2*i].z(),
+ gf->stl_curvatures[2*i+1].x(),gf->stl_curvatures[2*i+1].y(),gf->stl_curvatures[2*i+1].z());
+
}
for (size_t i=0;i<gf->stl_triangles.size()/3;i++){
fprintf(fp,"%d %d %d\n",gf->stl_triangles[3 * i + 0],gf->stl_triangles[3 * i + 1],gf->stl_triangles[3 * i + 2]);
diff --git a/Geo/GModelParametrize.cpp b/Geo/GModelParametrize.cpp
index 489941ebf44a87d64e0203aba1b2d4b62db0eca6..bbd18693916c155c76da270586591391f18b8123 100644
--- a/Geo/GModelParametrize.cpp
+++ b/Geo/GModelParametrize.cpp
@@ -79,6 +79,36 @@ static HXTStatus gmsh2hxt(GFace *gf, HXTMesh **pm,
return HXT_STATUS_OK;
}
+HXTStatus computeDiscreteCurvatures(GModel *gm, std::map<MVertex* , std::pair<SVector3, SVector3> > &C ){
+ C.clear();
+
+ for(GModel::fiter it = gm->firstFace(); it != gm->lastFace(); ++it) {
+ HXTMesh *m;
+ HXTEdges *edges;
+ double *nodalCurvatures;
+ double *crossField;
+ std::map<MVertex *, int> v2c;
+ std::vector<MVertex *> c2v;
+ gmsh2hxt(*it, &m, v2c, c2v);
+ HXT_CHECK(hxtEdgesCreate(m, &edges));
+ HXT_CHECK(hxtCurvatureRusinkiewicz(m, &nodalCurvatures, &crossField, edges, false));
+ for (size_t i=0;i<m->vertices.num;i++){
+ MVertex *v = c2v[i];
+ double *c = &nodalCurvatures[6*i];
+ SVector3 cMax (c[0],c[1],c[2]);
+ SVector3 cMin (c[3],c[4],c[5]);
+ std::pair<SVector3, SVector3> out = std::make_pair(cMax,cMin);
+ C.insert(std::make_pair(v,out));
+ }
+ HXT_CHECK(hxtEdgesDelete(&edges));
+ HXT_CHECK(hxtFree(&nodalCurvatures));
+ HXT_CHECK(hxtFree(&crossField));
+ HXT_CHECK(hxtMeshDelete(&m));
+ }
+ return HXT_STATUS_OK;
+}
+
+
void parametrizeAllGEdge(GModel *gm)
{
for(GModel::eiter it = gm->firstEdge(); it != gm->lastEdge(); ++it) {
@@ -87,17 +117,17 @@ void parametrizeAllGEdge(GModel *gm)
}
}
-int parametrizeAllGFace(GModel *gm)
+int parametrizeAllGFace(GModel *gm, std::map<MVertex* , std::pair<SVector3, SVector3> > *C)
{
int result = 0;
for(GModel::fiter it = gm->firstFace(); it != gm->lastFace(); ++it) {
discreteFace *df = dynamic_cast<discreteFace *>(*it);
- if(df) result += parametrizeGFace(df);
+ if(df) result += parametrizeGFace(df,C);
}
return result;
}
-int parametrizeGFace(discreteFace *gf)
+int parametrizeGFace(discreteFace *gf, std::map<MVertex* , std::pair<SVector3, SVector3> > *C)
{
int n = 1;
if(gf->triangles.empty()) return 0;
@@ -110,6 +140,7 @@ int parametrizeGFace(discreteFace *gf)
gmsh2hxt(gf, &m, v2c, c2v);
// double *nodalCurvatures;
// double *crossField;
+ // printf("A %d triangles\n",m->triangles.num);
HXT_CHECK(hxtEdgesCreate(m, &edges));
HXT_CHECK(hxtMeanValuesCreate(edges, ¶m));
HXT_CHECK(hxtMeanValuesCompute(param));
@@ -120,10 +151,20 @@ int parametrizeGFace(discreteFace *gf)
gf->stl_vertices_uv.resize(nv);
gf->stl_vertices_xyz.clear();
gf->stl_vertices_xyz.resize(nv);
+ gf->stl_curvatures.clear();
+ if (C) gf->stl_curvatures.resize(2*nv);
gf->stl_normals.clear();
gf->stl_normals.resize(nv);
-
+
for(int iv = 0; iv < nv; iv++) {
+ if (C){
+ MVertex *v = c2v[iv];
+ std::map<MVertex* , std::pair<SVector3, SVector3> >::iterator it = C->find(v);
+ if (it == C->end())Msg::Error("POINT %d NOT FOUND FOR COMPUTING CURVATURES",v->getNum());
+ // printf("%g %g %g\n",it->second.first.x(),it->second.first.y(),it->second.first.z());
+ gf->stl_curvatures[2*iv] = it->second.first;
+ gf->stl_curvatures[2*iv+1] = it->second.second;
+ }
gf->stl_vertices_uv[iv] = SPoint2(uvc[2*iv],uvc[2*iv+1]);
gf->stl_vertices_xyz[iv] =
@@ -138,14 +179,16 @@ int parametrizeGFace(discreteFace *gf)
gf->stl_triangles[3 * ie + 1] = m->triangles.node[3 * ie + 1];
gf->stl_triangles[3 * ie + 2] = m->triangles.node[3 * ie + 2];
}
+ // printf("A\n");
gf->fillVertexArray(false);
HXT_CHECK(hxtMeshDelete(&m));
HXT_CHECK(hxtEdgesDelete(&edges));
HXT_CHECK(hxtFree(&uvc));
+ // printf("B\n");
return 0;
}
#else
-int parametrizeGFace(GFace *gf)
+int parametrizeGFace(GFace *gf, std::map<std::pair<MVertex*,GFace*> , std::pair<SVector3, SVector3> > *C)
{
Msg::Error("Gmsh should be compiled against HXT for being able to compute "
"discrete parametrizations");
@@ -174,12 +217,27 @@ int isTriangulationParametrizable(const std::vector<MTriangle *> &t, int Nmax,
}
std::map<MEdge, int, Less_Edge>::iterator it = e.begin();
std::vector<MEdge> _bnd;
- for(; it != e.end(); ++it)
+ for(; it != e.end(); ++it){
if(it->second == 1) _bnd.push_back(it->first);
+ }
+
std::vector<std::vector<MVertex *> > vs;
- // printf("%d triangles %d edges %d vertices %d edges on
- // boundaries\n",t.size(),v.size(),e.size(),_bnd.size());
- if(!SortEdgeConsecutive(_bnd, vs)) return 2;
+ if(!SortEdgeConsecutive(_bnd, vs)) {
+ // we have vertices adjacent to more than 2 model edges
+ // FILE *f = fopen("bug.pos","w");
+ // fprintf(f,"View\"\"{\n");
+ // for (size_t i=0;i<t.size();i++){
+ // fprintf(f,"ST(%g,%g,%g,%g,%g,%g,%g,%g,%g){%d,%d,%d};\n",
+ // t[i]->getVertex(0)->x(),t[i]->getVertex(0)->y(),t[i]->getVertex(0)->z(),
+ // t[i]->getVertex(1)->x(),t[i]->getVertex(1)->y(),t[i]->getVertex(1)->z(),
+ // t[i]->getVertex(2)->x(),t[i]->getVertex(2)->y(),t[i]->getVertex(2)->z(),
+ // t[i]->getVertex(0)->getNum(),t[i]->getVertex(1)->getNum(),t[i]->getVertex(2)->getNum());
+ // }
+ // fprintf(f,"};\n");
+ // fclose(f);
+ // getchar();
+ return 2;
+ }
double lmax = 0;
for(size_t i = 0; i < vs.size(); i++) {
double li = 0;
@@ -189,7 +247,7 @@ int isTriangulationParametrizable(const std::vector<MTriangle *> &t, int Nmax,
lmax = std::max(li, lmax);
}
- if(ar * lmax * lmax < 2 * M_PI * surf) { return 2; }
+ if(ar * lmax * lmax < 2 * M_PI * surf) { return 4; }
int poincare =
t.size() - (2 * (v.size() - 1) - _bnd.size() + 2 * (vs.size() - 1));
@@ -206,6 +264,26 @@ void computeEdgeCut(GModel *gm, std::vector<MLine *> &cut,
int max_elems_per_cut)
{
GModel m;
+
+ // ----------------------------------------------------------------------------------
+ // STUPID FIX !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ for(GModel::fiter it = gm->firstFace(); it != gm->lastFace(); ++it) {
+ std::vector<MTriangle*> aa;
+ for(size_t i = 0; i < (*it)->triangles.size(); i++){
+ if ((*it)->triangles[i]->getVertex(0) == (*it)->triangles[i]->getVertex(1) ||
+ (*it)->triangles[i]->getVertex(0) == (*it)->triangles[i]->getVertex(2) ||
+ (*it)->triangles[i]->getVertex(1) == (*it)->triangles[i]->getVertex(2)){
+ // printf("TRIANGLE BAD %d %d %d\n",(*it)->triangles[i]->getVertex(0)->getNum()
+ // ,(*it)->triangles[i]->getVertex(1)->getNum()
+ // ,(*it)->triangles[i]->getVertex(2)->getNum());
+ }
+ else aa.push_back((*it)->triangles[i]);
+ }
+ (*it)->triangles = aa;
+ }
+ // END STUPID FIX !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ // ----------------------------------------------------------------------------------
+
for(GModel::fiter it = gm->firstFace(); it != gm->lastFace(); ++it) {
int part = 0;
if((*it)->triangles.empty()) continue;
diff --git a/Geo/GModelParametrize.h b/Geo/GModelParametrize.h
index 8cf198bd4ee7dc0a0c15bf64b463c5c1a5d2ee1b..22b93922d80fb8043222ff7b60f7046cf226a7c2 100644
--- a/Geo/GModelParametrize.h
+++ b/Geo/GModelParametrize.h
@@ -13,8 +13,9 @@ void computeEdgeCut(GModel *gm, std::vector<MLine *> &cut,
int max_elems_per_cut);
void computeNonManifoldEdges(GModel *gm, std::vector<MLine *> &cut,
bool addBoundary);
-int parametrizeAllGFace(GModel *gm);
+int parametrizeAllGFace(GModel *gm, std::map<MVertex* , std::pair<SVector3, SVector3> > *C=NULL);
void parametrizeAllGEdge(GModel *gm);
-int parametrizeGFace(discreteFace *gf);
+int parametrizeGFace(discreteFace *gf, std::map<MVertex* , std::pair<SVector3, SVector3> > *C=NULL);
+HXTStatus computeDiscreteCurvatures(GModel *gm, std::map<MVertex* , std::pair<SVector3, SVector3> > &C );
#endif
diff --git a/Geo/MEdge.cpp b/Geo/MEdge.cpp
index cbde5cfc662f7ad6c8c4a8cb24cb3e992f0d0d07..ff2f9b6a5ae640abb9727ef8903fb411e07ac46a 100644
--- a/Geo/MEdge.cpp
+++ b/Geo/MEdge.cpp
@@ -83,7 +83,7 @@ bool SortEdgeConsecutive(const std::vector<MEdge> &e,
else {
if(it0->second.second == NULL) { it0->second.second = v1; }
else {
- Msg::Error(
+ Msg::Debug(
"A list of edges is has points that are adjacet to 3 edges ! ");
return false;
}
@@ -93,7 +93,14 @@ bool SortEdgeConsecutive(const std::vector<MEdge> &e,
else {
if(it1->second.second == NULL) { it1->second.second = v0; }
else {
- Msg::Error("Wrong topology for a list of edges ");
+ Msg::Debug("Wrong topology for a list of edges ");
+ Msg::Debug("Vertex %d is adjacent to more than 2 vertices %d %d",v1->getNum(),it1->second.first->getNum(),it1->second.second->getNum());
+ // printf("Vertex %d is adjacent to more than 2 vertices %d %d\n",v1->getNum(),it1->second.first->getNum(),it1->second.second->getNum());
+ // for (size_t j=0;j<e.size();j++){
+ // printf("(%d %d)",e[j].getVertex(0)->getNum(),
+ // e[j].getVertex(1)->getNum());
+ // }
+ // printf("\n");
return false;
}
}
diff --git a/Geo/discreteEdge.cpp b/Geo/discreteEdge.cpp
index 982870c19f3c7bfceb53c6fffc53c873356c70b1..0c12436aad53f9ef943a181a353ce9dc1d3c9b7c 100644
--- a/Geo/discreteEdge.cpp
+++ b/Geo/discreteEdge.cpp
@@ -152,8 +152,18 @@ SVector3 discreteEdge::firstDer(double par) const
return der;
}
+SPoint2 discreteEdge::reparamOnFace(const GFace *face, double epar, int dir) const {
+ GPoint p = point(epar);
+ double guess[2];
+ GPoint ps = face->closestPoint(SPoint3(p.x(),p.y(),p.z()),guess);
+ // printf("garouclaboucazou %g %g %g %g %g %g\n",p.x(),p.y(),p.z(),ps.x(),ps.y(),ps.z());
+ return SPoint2(ps.u(),ps.v());
+}
+
+
double discreteEdge::curvature(double par) const
{
+ return 1.e-12;
double tLoc;
int iEdge;
if(_discretization.size() <= 3)
@@ -238,8 +248,14 @@ void discreteEdge::createGeometry()
if(lines.empty()) return;
if(!_discretization.empty()) return;
-
+
+ // printf("line %d\n",tag());
+ // for(std::size_t i = 0; i < lines.size(); i++) {
+ // printf("(%d %d)",lines[i]->getVertex(0)->getNum(),lines[i]->getVertex(1)->getNum());
+ // }
+ // printf("\n");
orderMLines();
+ // printf("done \n");
bool reverse = false;
if(getEndVertex())
diff --git a/Geo/discreteEdge.h b/Geo/discreteEdge.h
index 0365d161de04e2857e4bc0bb0c2dbaddc20df446..aa4d4c603284b3d6e281c1ad0fce5a9062fe0330 100644
--- a/Geo/discreteEdge.h
+++ b/Geo/discreteEdge.h
@@ -32,6 +32,7 @@ public:
virtual void mesh(bool verbose);
int minimumDrawSegments() const { return 2 * _pars.size(); }
virtual int minimumMeshSegments() const { return periodic(0) ? 3 : 2; }
+ virtual SPoint2 reparamOnFace(const GFace *face, double epar, int dir) const;
void setSplit(discreteEdge *e0, discreteEdge *e1, GVertex *vs)
{
_split[0] = e0;
diff --git a/Mesh/BackgroundMeshTools.cpp b/Mesh/BackgroundMeshTools.cpp
index 0a234f89671d85a7ddf4b0293b6077202c0e18b1..db69a2de1405a59114146f6e35fa8c4a32350d3a 100644
--- a/Mesh/BackgroundMeshTools.cpp
+++ b/Mesh/BackgroundMeshTools.cpp
@@ -14,12 +14,12 @@
static double max_surf_curvature(const GEdge *ge, double u)
{
- double val = 0;
+ double val = 0;
std::vector<GFace *> faces = ge->faces();
std::vector<GFace *>::iterator it = faces.begin();
while(it != faces.end()) {
SPoint2 par = ge->reparamOnFace((*it), u, 1);
- double cc = (*it)->curvature(par);
+ double cc = (*it)->curvatureMax(par);
val = std::max(cc, val);
++it;
}
diff --git a/Mesh/Generator.cpp b/Mesh/Generator.cpp
index 2292bbefdaf0a9455937c69ea238158cd529447d..c03e5d68d31fe84fe808aa7191899b03193dee88 100644
--- a/Mesh/Generator.cpp
+++ b/Mesh/Generator.cpp
@@ -482,10 +482,10 @@ static void Mesh2D(GModel *m)
for(GModel::fiter it = m->firstFace(); it != m->lastFace(); ++it) {
// STL remeshing is not yet thread-safe
- if((*it)->geomType() == GEntity::DiscreteSurface){
- if(static_cast<discreteFace *>(*it)->haveParametrization())
- Msg::SetNumThreads(1);
- }
+ // if((*it)->geomType() == GEntity::DiscreteSurface){
+ // if(static_cast<discreteFace *>(*it)->haveParametrization())
+ // Msg::SetNumThreads(1);
+ // }
// Frontal-Delaunay for quads and co are not yet thread-safe
if((*it)->getMeshingAlgo() == ALGO_2D_FRONTAL_QUAD ||
(*it)->getMeshingAlgo() == ALGO_2D_PACK_PRLGRMS ||
diff --git a/contrib/hxt/hxt_mean_values.c b/contrib/hxt/hxt_mean_values.c
index 1983923af9b9c713d0d5768ec3e2e758c900754a..adb2c4761018432884624790c8a097c581380c7b 100644
--- a/contrib/hxt/hxt_mean_values.c
+++ b/contrib/hxt/hxt_mean_values.c
@@ -34,7 +34,10 @@ HXTStatus hxtMeanValuesCreate(HXTEdges *edges, HXTMeanValues **meanValues)
// be moved into parametrization?
HXTBoundaries *boundaries;
+ printf("1 %d %d %d\n",mesh->triangles.node[0],mesh->triangles.node[1],mesh->triangles.node[2]);
+ printf("1 %d %d %d\n",mesh->triangles.node[3],mesh->triangles.node[4],mesh->triangles.node[5]);
HXT_CHECK(hxtEdgesSetBoundaries(edges, &boundaries));
+ printf("2\n");
map->boundaries = boundaries;
@@ -82,11 +85,12 @@ HXTStatus hxtMeanValuesCompute(HXTMeanValues *meanValues){
HXT_CHECK(hxtBoundariesGetNumberOfEdgesOfLineLoop(meanValues->boundaries,meanValues->hole[i],&c));
nby += c;
}
-
+
+#ifdef __FILLERGAP_
uint32_t *fillergap;
HXT_CHECK(hxtMalloc(&fillergap, 3*(nTriangles+nby)*sizeof(uint32_t)));
memcpy(fillergap, mesh->triangles.node, 3*nTriangles*sizeof(uint32_t));
-
+
int s=0;
for(int i=0; i<meanValues->nHoles; i++){
uint32_t *cedgs;
@@ -100,20 +104,27 @@ HXTStatus hxtMeanValuesCompute(HXTMeanValues *meanValues){
}
s += n_;
}
+#endif // __FILLERGAP_
+
/*
for(int i=0; i<nTriangles+nby; i++)
printf("tri %d (nby=%d) \t %u %u %u\n",i,nby,fillergap[3*i+0],fillergap[3*i+1],fillergap[3*i+2]);
*/
-
+
HXTLinearSystem *sys;
-
+
+#ifdef __FILLERGAP_
HXT_CHECK(hxtLinearSystemCreateLU(&sys,nTriangles+nby,3,1,fillergap));
-
+#else
+ HXT_CHECK(hxtLinearSystemCreateLU(&sys,nTriangles,3,1,mesh->triangles.node));
+#endif
+
+
uint32_t *flag;
HXT_CHECK(hxtMalloc(&flag,nNodes*sizeof(uint32_t)));
for(uint32_t ii=0; ii<nNodes; ii++)
flag[ii] = 0;
-
+
/* boundary conditions */
double *uv = meanValues->uv;
int n;
@@ -130,28 +141,28 @@ HXTStatus hxtMeanValuesCompute(HXTMeanValues *meanValues){
uv[2*edges->node[2*edges_ll[i]+0]+1] = sin(angle);
currentLength += hxtEdgesLength(meanValues->initialEdges, edges_ll[i]);
}
-
+
double *U, *V, *Urhs, *Vrhs;
//printf("allocation: %d\n",nNodes+meanValues->nHoles);
HXT_CHECK( hxtMalloc(&U,(nNodes+meanValues->nHoles)*sizeof(double)) );
HXT_CHECK( hxtMalloc(&V,(nNodes+meanValues->nHoles)*sizeof(double)) );
HXT_CHECK( hxtMalloc(&Urhs,(nNodes+meanValues->nHoles)*sizeof(double)) );
HXT_CHECK( hxtMalloc(&Vrhs,(nNodes+meanValues->nHoles)*sizeof(double)) );
-
-
+
+
// init linear system
HXT_CHECK(hxtLinearSystemZeroMatrix(sys));
for(uint32_t i=0; i<(nNodes+meanValues->nHoles); i++){
Urhs[i] = 0.;
Vrhs[i] = 0.;
}
-
+
// setting linear system
for(uint32_t ie=0; ie<edges->numEdges; ie++){
int ik[2] = {-1,-1};
uint64_t *tri = edges->edg2tri + 2*ie;
-
+
// gather local edge index
for(int it=0; it<2; it++){
if(tri[it]==(uint64_t)-1)
@@ -163,21 +174,21 @@ HXTStatus hxtMeanValuesCompute(HXTMeanValues *meanValues){
}
}
}
-
+
for(int ij=0; ij<2; ij++){
-
+
uint32_t v0 = edges->node[2*ie+ij];
uint32_t v1 = edges->node[2*ie+(1-ij)];
-
+
if (flag[v0]==1){//boundary nodes/conditons
- HXT_CHECK(hxtLinearSystemSetMatrixRowIdentity(sys,v0,0));
+ HXT_CHECK(hxtLinearSystemSetMatrixRowIdentity(sys,v0,0));
HXT_CHECK(hxtLinearSystemSetRhsEntry(sys,Urhs, v0,0, uv[2*v0+0]));
HXT_CHECK(hxtLinearSystemSetRhsEntry(sys,Vrhs, v0,0, uv[2*v0+1]));
}
else {// inner node
- double e[3] = {mesh->vertices.coord[4*v1+0]-mesh->vertices.coord[4*v0+0],mesh->vertices.coord[4*v1+1]-mesh->vertices.coord[4*v0+1],mesh->vertices.coord[4*v1+2]-mesh->vertices.coord[4*v0+2]};
- double ne = sqrt(e[0]*e[0]+e[1]*e[1]+e[2]*e[2]);
-
+ double e[3] = {mesh->vertices.coord[4*v1+0]-mesh->vertices.coord[4*v0+0],mesh->vertices.coord[4*v1+1]-mesh->vertices.coord[4*v0+1],mesh->vertices.coord[4*v1+2]-mesh->vertices.coord[4*v0+2]};
+ double ne = sqrt(e[0]*e[0]+e[1]*e[1]+e[2]*e[2]);
+
uint32_t vLeft = mesh->triangles.node[3*tri[0] + (ik[0]+2)%3];
double a[3] = {mesh->vertices.coord[4*vLeft+0]-mesh->vertices.coord[4*v0+0],mesh->vertices.coord[4*vLeft+1]-mesh->vertices.coord[4*v0+1],mesh->vertices.coord[4*vLeft+2]-mesh->vertices.coord[4*v0+2]};
double na = sqrt(a[0]*a[0]+a[1]*a[1]+a[2]*a[2]);
@@ -199,9 +210,9 @@ HXTStatus hxtMeanValuesCompute(HXTMeanValues *meanValues){
}
}
}// end for int ie
+
-
-
+#ifdef __FILLERGAP_
// holes are filled somehow
for(int i=0; i<meanValues->nHoles; i++){
@@ -258,6 +269,8 @@ HXTStatus hxtMeanValuesCompute(HXTMeanValues *meanValues){
HXT_CHECK(hxtFree(&c));
HXT_CHECK(hxtFree(&d));
}
+#endif
+
HXT_CHECK(hxtLinearSystemSolve(sys,Urhs,U));
HXT_CHECK(hxtLinearSystemSolve(sys,Vrhs,V));
for(uint32_t i=0; i<nNodes; i++){
@@ -265,7 +278,9 @@ HXTStatus hxtMeanValuesCompute(HXTMeanValues *meanValues){
uv[2*i+1] = V[i] ;
}
+#ifdef __FILLERGAP_
HXT_CHECK(hxtFree(&fillergap));
+#endif
HXT_CHECK(hxtFree(&flag));
HXT_CHECK(hxtFree(&U));
HXT_CHECK(hxtFree(&V));