diff --git a/Mesh/CenterlineField.cpp b/Mesh/CenterlineField.cpp
index c33087796f1fed3377f4156c5d28b65159f69006..03899fca412344a61247b25b431c05218d538a6e 100644
--- a/Mesh/CenterlineField.cpp
+++ b/Mesh/CenterlineField.cpp
@@ -11,6 +11,8 @@
#include <vector>
#include <map>
+#include <set>
+#include <list>
#include <string>
#include <fstream>
#include <sstream>
@@ -25,6 +27,129 @@
#if defined(HAVE_ANN)
#include <ANN/ANN.h>
+#include <algorithm>
+
+void erase(std::vector<MLine*>& lines, MLine* l) {
+ std::vector<MLine*>::iterator it = std::remove(lines.begin(), lines.end(), l);
+ lines.erase(it, lines.end());
+}
+double computeLength(std::vector<MLine*> lines){
+
+ double length= 0.0;
+ for (int i = 0; i< lines.size(); i++){
+ length += lines[i]->getInnerRadius();
+ }
+ return length;
+}
+
+void orderMLines(std::vector<MLine*> &lines) { //, std::set<MVertex*> & junctions, MVertex *vB, MVertex *vE){
+
+ std::vector<MLine*> _m;
+ std::list<MLine*> segments;
+ std::map<MVertex*, MLine*> boundv;
+ std::vector<int> _orientation;
+
+ // store all lines in a list : segments
+ for (unsigned int i = 0; i < lines.size(); i++){
+ segments.push_back(lines[i]);
+ }
+
+ // find a lonely MLine
+ for (std::list<MLine*>::iterator it = segments.begin();
+ it != segments.end(); ++it){
+ MVertex *vL = (*it)->getVertex(0);
+ MVertex *vR = (*it)->getVertex(1);
+ std::map<MVertex*,MLine*>::iterator it1 = boundv.find(vL);
+ if (it1 == boundv.end()) boundv.insert(std::make_pair(vL,*it));
+ else boundv.erase(it1);
+ std::map<MVertex*,MLine*>::iterator it2 = boundv.find(vR);
+ if (it2 == boundv.end()) boundv.insert(std::make_pair(vR,*it));
+ else boundv.erase(it2);
+ }
+
+ // find the first MLine and erase it from the list segments
+ MLine *firstLine;
+ if (boundv.size() == 2){ // non periodic
+ firstLine = (boundv.begin())->second;
+ for (std::list<MLine*>::iterator it = segments.begin();
+ it != segments.end(); ++it){
+ if (*it == firstLine){
+ segments.erase(it);
+ break;
+ }
+ }
+ }
+ else if (boundv.size() == 0){ // periodic
+ firstLine = *(segments.begin());
+ segments.erase(segments.begin());
+ }
+ else{
+ Msg::Error("line is wrong (it has %d end points)",
+ boundv.size());
+ }
+
+ // loop over all segments to order segments and store it in the list _m
+ _m.push_back(firstLine);
+ _orientation.push_back(1);
+ MVertex *first = _m[0]->getVertex(0);
+ MVertex *last = _m[0]->getVertex(1);
+ while (first != last){
+ if (segments.empty())break;
+ bool found = false;
+ for (std::list<MLine*>::iterator it = segments.begin();
+ it != segments.end(); ++it){
+ MLine *e = *it;
+ if (e->getVertex(0) == last){
+ _m.push_back(e);
+ segments.erase(it);
+ _orientation.push_back(1);
+ last = e->getVertex(1);
+ found = true;
+ break;
+ }
+ else if (e->getVertex(1) == last){
+ _m.push_back(e);
+ segments.erase(it);
+ _orientation.push_back(0);
+ last = e->getVertex(0);
+ found = true;
+ break;
+ }
+ }
+ if (!found && _orientation[0]==1){ //reverse orientation of first Line
+ if (_m.size() == 1 ){
+ MVertex *temp = first;
+ first = last;
+ last = temp;
+ _orientation[0] = 0;
+ }
+ else {
+ Msg::Error("lines is wrong");
+ return;
+ }
+ }
+ }
+
+ //lines is now a list of ordered MLines
+ lines = _m;
+
+ //special case reverse orientation
+ if (lines.size() < 2) return;
+ if (_orientation[0] && lines[0]->getVertex(1) != lines[1]->getVertex(1)
+ && lines[0]->getVertex(1) != lines[1]->getVertex(0)){
+ for (unsigned int i = 0; i < lines.size(); i++)
+ _orientation[i] = !_orientation[i];
+ }
+
+ // if (junctions.find(lines[0]->getVertex(0)) != junctions.end()) vB = lines[0]->getVertex(0);
+ // else if (junctions.find(lines[0]->getVertex(1)) != junctions.end()) vB = lines[0]->getVertex(1);
+ // else printf("no vB junc found %d %d \n", lines[0]->getVertex(0)->getNum(), lines[0]->getVertex(1)->getNum());
+ // if (junctions.find(lines[lines.size()-1]->getVertex(0)) != junctions.end()) vE = lines[lines.size()-1]->getVertex(0);
+ // else if (junctions.find(lines[lines.size()-1]->getVertex(1)) != junctions.end()) vE = lines[lines.size()-1]->getVertex(1);
+ // else printf("no vE junc found %d %d \n", lines[0]->getVertex(0)->getNum(), lines[0]->getVertex(1)->getNum());
+ // printf("in order vB= %d =%d \n", vB->getNum(), vE->getNum());
+}
+
Centerline::Centerline(std::string fileName): kdtree(0), nodes(0){
index = new ANNidx[1];
@@ -34,7 +159,6 @@ Centerline::Centerline(std::string fileName): kdtree(0), nodes(0){
importFile(fileName);
buildKdTree();
-
}
Centerline::Centerline(): kdtree(0), nodes(0){
@@ -64,44 +188,158 @@ void Centerline::importFile(std::string fileName){
std::vector<GEntity*> entities ;
mod->getEntities(entities) ;
-
+
+ int maxN = 0.0;
for(unsigned int i = 0; i < entities.size(); i++){
if( entities[i]->dim() == 1){
for(unsigned int ele = 0; ele < entities[i]->getNumMeshElements(); ele++){
MLine *l = (MLine*) entities[i]->getMeshElement(ele);
MVertex *v0 = l->getVertex(0);
MVertex *v1 = l->getVertex(1);
- std::map<MVertex*, int>::iterator it0 = color.find(v0);
- std::map<MVertex*, int>::iterator it1 = color.find(v1);
- if (it0 == color.end() || it1 == color.end()) lines.push_back(l);
- if (it0 == color.end()) color.insert(std::make_pair(v0, entities[i]->tag()));
- if (it1 == color.end()) color.insert(std::make_pair(v1, entities[i]->tag()));
+ std::map<MVertex*, int>::iterator it0 = colorp.find(v0);
+ std::map<MVertex*, int>::iterator it1 = colorp.find(v1);
+ if (it0 == colorp.end() || it1 == colorp.end()){
+ lines.push_back(l);
+ colorl.insert(std::make_pair(l, entities[i]->tag()));
+ maxN = std::max(maxN, entities[i]->tag());
+ }
+ if (it0 == colorp.end()) colorp.insert(std::make_pair(v0, entities[i]->tag()));
+ if (it1 == colorp.end()) colorp.insert(std::make_pair(v1, entities[i]->tag()));
}
}
}
+
+ //splitEdges
+ splitEdges(maxN);
+
current->setAsCurrent();
}
+
+void Centerline::splitEdges(int maxN){
+
+ //sort colored lines and create edges
+ std::vector<std::vector<MLine*> > color_edges;
+ color_edges.resize(maxN);
+ std::multiset<MVertex*> allV;
+ std::map<MLine*, int>::iterator itl = colorl.begin();
+ while (itl != colorl.end()){
+ int color = itl->second;
+ MLine* l = itl->first;
+ allV.insert(l->getVertex(0));
+ allV.insert(l->getVertex(1));
+ color_edges[color-1].push_back(l);
+ itl++;
+ }
+
+ //detect junctions
+ std::multiset<MVertex*>::iterator it = allV.begin();
+ for ( ; it != allV.end(); ++it){
+ if (allV.count(*it) != 2) {
+ junctions.insert(*it);
+ }
+ }
+
+ //split edges
+ int tag = 0;
+ for(unsigned int i = 0; i < color_edges.size(); ++i){
+ std::vector<MLine*> lines = color_edges[i];
+ while (!lines.empty()) {
+ std::vector<MLine*> myLines;
+ std::vector<MLine*>::iterator itl = lines.begin();
+ MVertex *vB = (*itl)->getVertex(0);
+ MVertex *vE = (*itl)->getVertex(1);
+ myLines.push_back(*itl);
+ erase(lines, *itl);
+ itl = lines.begin();
+ while ( !( junctions.find(vE) != junctions.end() &&
+ junctions.find(vB) != junctions.end()) ) {
+ MVertex *v1 = (*itl)->getVertex(0);
+ MVertex *v2 = (*itl)->getVertex(1);
+ if (v1 == vE){
+ myLines.push_back(*itl);
+ erase(lines, *itl);
+ itl = lines.begin();
+ vE = v2;
+ }
+ else if ( v2 == vE){
+ myLines.push_back(*itl);
+ erase(lines, *itl);
+ itl = lines.begin();
+ vE = v1;
+ }
+ else if ( v1 == vB){
+ myLines.push_back(*itl);
+ erase(lines, *itl);
+ itl = lines.begin();
+ vB = v2;
+ }
+ else if ( v2 == vB){
+ myLines.push_back(*itl);
+ erase(lines, *itl);
+ itl = lines.begin();
+ vB = v1;
+ }
+ else itl++;
+ }
+ if (vB == vE) { Msg::Error("Begin and end points branch are the same \n");break;}
+ orderMLines(myLines); //, junctions, vBB, vEE);
+ std::vector<Branch> children;
+ Branch newBranch ={ tag++, myLines, computeLength(myLines), vB, vE, children};
+ //printf("VB = %d %d \n", vB->getNum(), vE->getNum());
+ edges.push_back(newBranch) ;
+ }
+ }
+
+ //create new edges
+ printf("edges =%d new =%d \n", color_edges.size(), edges.size());
+
+ //create children
+ for(unsigned int i = 0; i < edges.size(); ++i) {
+ MVertex *vE = edges[i].vE;
+ std::vector<Branch> myChildren;
+ for (std::vector<Branch>::iterator it = edges.begin(); it != edges.end(); ++it){
+ Branch myBranch = *it;
+ if (myBranch.vB == vE) myChildren.push_back(myBranch);
+ }
+ edges[i].children = myChildren;
+ }
+
+ //print info
+ // for(unsigned int i = 0; i < edges.size(); ++i) {
+ // printf("EDGE =%d tag=%d length = %g childs = %d \n", i, edges[i].tag, edges[i].length, edges[i].children.size());
+ // for(unsigned int j = 0; j < edges[i].children.size(); ++j) {
+ // printf("children (%d) =%d \n", edges[i].children.size(), edges[i].children[j].tag);
+ // }
+ // }
+
+ // std::set<MVertex*>::iterator itj = junctions.begin();
+ // for ( ; itj != junctions.end(); ++itj){
+ // MVertex *v = *itj;
+ // printf("JUNC = %d \n", v->getNum());
+ // }
+
+}
+
void Centerline::buildKdTree(){
FILE * f = fopen("myPOINTS.pos","w");
fprintf(f, "View \"\"{\n");
- int nbPL = 10;
+ int nbPL = 3; //10 points per line
int nbNodes = (lines.size()+1) + (nbPL*lines.size());
nodes = annAllocPts(nbNodes, 3);
int ind = 0;
- std::map<MVertex*,int >::iterator itmap = color.begin();
- while (itmap != color.end()){
- MVertex *v = itmap->first;
+ std::map<MVertex*, int>::iterator itp = colorp.begin();
+ while (itp != colorp.end()){
+ MVertex *v = itp->first;
nodes[ind][0] = v->x();
nodes[ind][1] = v->y();
nodes[ind][2] = v->z();
- itmap++; ind++;
+ itp++; ind++;
}
- //10 points per line
- for(unsigned int k = 0; k < lines.size(); ++k){
+ for(unsigned int k = 0; k < lines.size(); ++k){
MVertex *v0 = lines[k]->getVertex(0);
MVertex *v1 = lines[k]->getVertex(1);
SVector3 P0(v0->x(),v0->y(), v0->z());
@@ -155,22 +393,70 @@ void Centerline::operator() (double x, double y, double z, SMetric3 &metr, GEnt
void Centerline::printSplit() const{
- FILE * f = fopen("centerlineSPLIT.pos","w");
+ FILE * f = fopen("mySPLIT.pos","w");
fprintf(f, "View \"\"{\n");
- for(unsigned int i = 0; i < lines.size(); ++i){
- MLine *l = lines[i];
- std::map<MVertex*,int>::const_iterator it0 =
- color.find(l->getVertex(0));
- std::map<MVertex*,int>::const_iterator it1 =
- color.find(l->getVertex(1));
+ // for(unsigned int i = 0; i < lines.size(); ++i){
+ // MLine *l = lines[i];
+ // std::map<MLine*,int>::const_iterator itc = colorl.find(l);
+ // std::map<MVertex*,int>::const_iterator it0 =
+ // colorp.find(l->getVertex(0));
+ // std::map<MVertex*,int>::const_iterator it1 =
+ // colorp.find(l->getVertex(1));
+ // fprintf(f, "SL(%g,%g,%g,%g,%g,%g){%g,%g};\n",
+ // l->getVertex(0)->x(), l->getVertex(0)->y(), l->getVertex(0)->z(),
+ // l->getVertex(1)->x(), l->getVertex(1)->y(), l->getVertex(1)->z(),
+ // (double)itc->second,(double)itc->second);
+ // }
+ // std::map<MVertex*,int>::const_iterator itp = colorp.begin();
+ // while (itp != colorp.end()){
+ // MVertex *v = itp->first;
+ // fprintf(f, "SP(%g,%g,%g){%g};\n",
+ // v->x(), v->y(), v->z(),
+ // 0.0);
+ // itp++;
+ // }
+
+ for(unsigned int i = 0; i < edges.size(); ++i){
+ std::vector<MLine*> lines = edges[i].lines;
+ for(unsigned int k = 0; k < lines.size(); ++k){
+ MLine *l = lines[k];
fprintf(f, "SL(%g,%g,%g,%g,%g,%g){%g,%g};\n",
- l->getVertex(0)->x(), l->getVertex(0)->y(), l->getVertex(0)->z(),
- l->getVertex(1)->x(), l->getVertex(1)->y(), l->getVertex(1)->z(),
- (double)it0->second,(double)it1->second);
- }
+ l->getVertex(0)->x(), l->getVertex(0)->y(), l->getVertex(0)->z(),
+ l->getVertex(1)->x(), l->getVertex(1)->y(), l->getVertex(1)->z(),
+ (double)i, (double)i);
+ }
+ }
+
fprintf(f,"};\n");
fclose(f);
+
+ FILE * f2 = fopen("myCOLORS.pos","w");
+ fprintf(f2, "View \"\"{\n");
+ std::map<MVertex*, int>::const_iterator itp = colorp.begin();
+ while (itp != colorp.end()){
+ MVertex *v = itp->first;
+ fprintf(f2, "SP(%g,%g,%g){%g};\n",
+ v->x(), v->y(), v->z(),
+ (double)v->getNum());
+ itp++;
+ }
+ fprintf(f2,"};\n");
+ fclose(f2);
+
+ FILE * f3 = fopen("myJUNCTIONS.pos","w");
+ fprintf(f3, "View \"\"{\n");
+ std::set<MVertex*>::const_iterator itj = junctions.begin();
+ while (itj != junctions.end()){
+ MVertex *v = *itj;
+ fprintf(f3, "SP(%g,%g,%g){%g};\n",
+ v->x(), v->y(), v->z(),
+ (double)v->getNum());
+ itj++;
+ }
+ fprintf(f3,"};\n");
+ fclose(f3);
+
}
diff --git a/Mesh/CenterlineField.h b/Mesh/CenterlineField.h
index ceb2c4ddfd4f8f8f484ec1894125e7ea51ea1f17..79f60fa1d96394274f0fae92391a674feafee864 100644
--- a/Mesh/CenterlineField.h
+++ b/Mesh/CenterlineField.h
@@ -11,6 +11,7 @@
#include <vector>
#include <map>
+#include <set>
#include <string>
#include "Field.h"
class GModel;
@@ -22,6 +23,15 @@ class GEntity;
#include <ANN/ANN.h>
class ANNkd_tree;
+struct Branch{
+ int tag;
+ std::vector<MLine*> lines;
+ double length;
+ MVertex *vB;
+ MVertex *vE;
+ std::vector<Branch> children;
+};
+
class Centerline : public Field{
ANNkd_tree *kdtree;
@@ -37,6 +47,7 @@ class Centerline : public Field{
~Centerline();
void importFile(std::string fileName);
+ void splitEdges(int maxN);
virtual bool isotropic () const {return false;}
virtual const char *getName()
@@ -58,7 +69,11 @@ class Centerline : public Field{
void buildKdTree();
std::vector<MLine*> lines;
- std::map<MVertex*,int> color;
+ std::vector<Branch> edges;
+ std::map<MVertex*,int> colorp;
+ std::map<MLine*,int> colorl;
+
+ std::set<MVertex*> junctions;
};
#endif