diff --git a/Geo/discreteDiskFace.cpp b/Geo/discreteDiskFace.cpp
index 017b76c684ffc53b558b379d32198ffcca3695e7..d8a00a4727f7897e726afcfdcbd3a1ce141df8bc 100644
--- a/Geo/discreteDiskFace.cpp
+++ b/Geo/discreteDiskFace.cpp
@@ -203,13 +203,14 @@ static bool orderVertices(const double &tot_length, const std::vector<MVertex*>
(next->z() - first->z()) * (next->z() - first->z()) );
coord.push_back(coord[coord.size()-1] + length / tot_length);
-
+
first = next;
}
return true;
}
+/*BUILDER*/
discreteDiskFace::discreteDiskFace(GFace *gf, triangulation* diskTriangulation,
int p, std::vector<GFace*> *CAD) :
GFace(gf->model(),123), _parent (gf), _ddft(NULL), oct(NULL)
@@ -280,9 +281,8 @@ discreteDiskFace::discreteDiskFace(GFace *gf, triangulation* diskTriangulation,
geoTriangulation = new triangulation(discrete_triangles,gf);
allNodes = geoTriangulation->vert;
- _U0 = geoTriangulation->bord.rbegin()->second;
_totLength = geoTriangulation->bord.rbegin()->first;
-
+ _U0 = geoTriangulation->bord.rbegin()->second;
orderVertices(_totLength, _U0, _coords);
@@ -290,14 +290,14 @@ discreteDiskFace::discreteDiskFace(GFace *gf, triangulation* diskTriangulation,
buildOct(CAD);
if (!checkOrientationUV()){
- Msg::Info("discreteDIskFace:: parametrization is not one-to-one; fixing "
+ Msg::Info("discreteDiskFace:: parametrization is not one-to-one; fixing "
"the discrete system.");
parametrize(true);
buildOct(CAD);
}
-
putOnView();
}
+/*end BUILDER*/
discreteDiskFace::~discreteDiskFace()
{
@@ -356,8 +356,14 @@ bool discreteDiskFace::parametrize(bool one2one) const
for(size_t i = 0; i < _U0.size(); i++){
MVertex *v = _U0[i];
const double theta = 2 * M_PI * _coords[i];
- myAssemblerU.fixVertex(v, 0, 1, cos(theta));
- myAssemblerV.fixVertex(v, 0, 1, sin(theta));
+ if(i%_order==0){
+ myAssemblerU.fixVertex(v, 0, 1, cos(theta));
+ myAssemblerV.fixVertex(v, 0, 1, sin(theta));
+ }
+ else{//#TEST
+ myAssemblerU.fixVertex(v, 0, 1,1./_order*((_order-(i%_order))*cos(2*M_PI*_coords[i-(i%_order)])+(i%_order)*cos(2*M_PI*_coords[i+(_order-(i%_order))])));
+ myAssemblerV.fixVertex(v, 0, 1,1./_order*((_order-(i%_order))*sin(2*M_PI*_coords[i-(i%_order)])+(i%_order)*sin(2*M_PI*_coords[i+(_order-(i%_order))])));
+ }
}
for(size_t i = 0; i < discrete_triangles.size(); ++i){
diff --git a/Geo/discreteDiskFace.h b/Geo/discreteDiskFace.h
index 20124e03a95015b4da7c368f7f700fd011bdd48c..bf333efa787ba241f98ae5324964e86a529add39 100644
--- a/Geo/discreteDiskFace.h
+++ b/Geo/discreteDiskFace.h
@@ -57,8 +57,10 @@ class triangulation {
GFace *gf;
int idNum; // number of identification, for hashing purposes
- int genus()
- {
+ std::list<GEdge*> my_GEdges;
+
+ //methods
+ int genus(){
return ( -vert.size() + ed2tri.size() - tri.size() + 2 - bord.size() )/2;
}
diff --git a/Geo/discreteFace.cpp b/Geo/discreteFace.cpp
index e9a69c260501d39828fb1dca0c664e064268ca04..3edbdf2455f24e92194fdb9fe27495d65e00e87c 100644
--- a/Geo/discreteFace.cpp
+++ b/Geo/discreteFace.cpp
@@ -105,12 +105,11 @@ void discreteFace::secondDer(const SPoint2 ¶m,
return;
}
-// FIXME PAB ----> really create an atlas !!!!!!!!!!!!!!!
void discreteFace::createGeometry()
{
checkAndFixOrientation();
- int order = 1;
+ int order = 2;
int nPart = 2;
std::vector<triangulation*> part;
part.resize(nPart);
@@ -122,7 +121,11 @@ void discreteFace::createGeometry()
std::stack<triangulation*> toSplit;
std::vector<triangulation*> toParam;
std::vector<MElement*> tem(triangles.begin(),triangles.end());
- toSplit.push(new triangulation(tem,this));
+
+ triangulation* init = new triangulation(tem,this);
+ init->my_GEdges = l_edges;
+
+ toSplit.push(init);
if((toSplit.top())->genus()!=0){
while( !toSplit.empty()){
@@ -139,6 +142,7 @@ void discreteFace::createGeometry()
else{
toParam.push_back(part[i]);
part[i]->idNum=id++;
+ // discreteEdges
}
}// end for i
}
@@ -148,10 +152,11 @@ void discreteFace::createGeometry()
toParam.push_back(toSplit.top());
toSplit.top()->idNum=id++;
}
-
- for(unsigned int i=0; i<toParam.size(); i++){
+
+ for(unsigned int i=0; i<toParam.size(); i++){
+ std::vector<MElement*> mytri = toParam[i]->tri;
discreteDiskFace *df = new discreteDiskFace (this,toParam[i], order,(_CAD.empty() ? NULL : &_CAD));
- df->replaceEdges(l_edges);// #FIXME
+ df->replaceEdges(toParam[i]->my_GEdges);
_atlas.push_back(df);
}
#endif
@@ -215,7 +220,7 @@ void discreteFace::checkAndFixOrientation(){
// first of all, all the triangles have to be oriented in the same way
std::map<MEdge,std::vector<MElement*>,Less_Edge> ed2tri; // edge to 1 or 2 triangle(s)
-
+
for(unsigned int i = 0; i < triangles.size(); ++i){
MElement *e = triangles[i];
for(int j = 0; j < e->getNumEdges() ; j++){
@@ -281,7 +286,7 @@ void discreteFace::checkAndFixOrientation(){
if (!(current->getVertex(k!=2 ?k+1 : 0 ) == neigs[i]->getVertex(j!=0 ? j-1 : 2) ||
current->getVertex(k!=0 ?k-1 : 2 ) == neigs[i]->getVertex(j!=2 ? j+1 : 0))){
neigs[i]->reverse();
- Msg::Info("discreteDiskFace: triangle %d has been reoriented.",neigs[i]->getNum());
+ Msg::Info("discreteFace: triangle %d has been reoriented.",neigs[i]->getNum());
}
break;
}
@@ -315,8 +320,8 @@ void discreteFace::split(triangulation* trian,std::vector<triangulation*> &parti
MEdge ed = current->getEdge(j);
int nEd = trian->ed2tri[ed].size();
-
- if (nEd > 1){
+
+ if (nEd > 1){
std::vector<int> local = trian->ed2tri[ed];
nbh[idx[i]+temp] = local[0] == i ? local[1] : local[0];
temp++;
@@ -333,17 +338,138 @@ void discreteFace::split(triangulation* trian,std::vector<triangulation*> &parti
std::vector<int> part;
part.resize(nVertex);
int zero = 0;
- METIS_PartGraphRecursive(&nVertex,&idx[0],&nbh[0],NULL,NULL,&zero,&zero,&nPartitions,&zero,&edgeCut,&part[0]);
- // CONNEC
+ METIS_PartGraphRecursive(&nVertex,&(idx[0]),&(nbh[0]),NULL,NULL,&zero,&zero,&nPartitions,&zero,&edgeCut,&(part[0]));
+
std::vector<std::vector<MElement*> > elem;
elem.resize(nPartitions);
for(int i=0; i<nVertex; i++)
elem[part[i]].push_back(trian->tri[i]);
-
+
for(int i=0; i<nPartitions; i++)
partition[i] = new triangulation(elem[i],this);
+
+
+ for(int i=0; i<nPartitions; i++){// setting GEdge's
+
+ //"update" old GEdge's
+ std::list<GEdge*> upGed = partition[i]->my_GEdges;// old GEdge's list
+ for(std::list<GEdge*>::iterator cGedg=trian->my_GEdges.begin(); cGedg!=trian->my_GEdges.end(); cGedg++){// GEdge by GEdge
+ std::set<MVertex*> first, last;
+ double lc = 0.;
+ std::vector<MLine*> upMli;
+ for(unsigned int k=0; k<(*cGedg)->lines.size(); k++){// MLine by MLine (for a GEdge)
+ MLine* cMli = (*cGedg)->lines[k];
+ std::set<MVertex*>::iterator mvit0 = partition[i]->vert.find(cMli->getVertex(0));
+ std::set<MVertex*>::iterator mvit1 = partition[i]->vert.find(cMli->getVertex(1));
+ if(mvit0 != partition[i]->vert.end() && mvit1 != partition[i]->vert.end()){// check if the MLine belong to this part
+
+ if( first.find(cMli->getVertex(1)) != first.end() )
+ first.erase(cMli->getVertex(1));
+ else
+ last.insert(cMli->getVertex(1));
+
+ if( last.find(cMli->getVertex(0)) != last.end() )
+ last.erase(cMli->getVertex(0));
+ else
+ first.insert(cMli->getVertex(0));
+
+ upMli.push_back(cMli);
+ lc += cMli->getLength();
+
+ }
+ }// end for k
+ GVertex *v0, *v1;
+ if(first.empty()){
+ v0 = this->model()->addVertex(upMli[0]->getVertex(0)->x(),upMli[0]->getVertex(0)->y(),upMli[0]->getVertex(0)->z(),1.);//#fixme
+ v1 = v0;
+ }
+ else{
+ v0 = this->model()->addVertex((*(first.begin()))->x(),(*(first.begin()))->y(),(*(first.begin()))->z(),1.);//#fixme
+ v1 = this->model()->addVertex((*(last.begin()))->x(),(*(last.begin()))->y(),(*(last.begin()))->z(),1.);
+ }
+ upGed.push_back(this->model()->addLine(v0,v1));
+ upGed.back()->lines = upMli;
+ }// end for j
+
+
+ // create new GEdge's
+ std::set<MEdge,Less_Edge> bordi = partition[i]->borderEdg;// edges defining the border(s) of the i-th new triangulation
+ for(int ii=i+1; ii<nPartitions; ii++){// compare to the ii-th ones, with ii > i since ii < i have already been compared
+ std::map<MVertex*,MLine*> v02edg;//first vertex of the MLine
+ std::set<MVertex*> first, last;
+ for(std::set<MEdge,Less_Edge>::iterator ie = bordi.begin();
+ ie != bordi.end(); ++ie){// MEdge by MEdge of the i-th triangulation border(s)
+
+ std::set<MEdge,Less_Edge> bii = partition[ii]->borderEdg;// edges defining the border(s) of the ii-th new triangulation
+ if(bii.find(*ie)!=bii.end()){// if the border edge is common to both triangulations . . .
+
+ v02edg[ie->getVertex(0)] = new MLine(ie->getVertex(0),ie->getVertex(1));// . . . a new MLine is created
+
+ // identifying first and last vertices of GEdge's, if any
+ if( first.find(ie->getVertex(1)) != first.end() )
+ first.erase(ie->getVertex(1));
+ else
+ last.insert(ie->getVertex(1));
+ if( last.find(ie->getVertex(0)) != last.end() )
+ last.erase(ie->getVertex(0));
+ else
+ first.insert(ie->getVertex(0));
+ }
+
+ }//end for ie
+
+
+ while(!v02edg.empty()){// creation of the GEdge's from new MLine's
+
+ GVertex *v0, *v1;
+
+ while(!first.empty()){// starting with nonloop GEdge's
+
+ std::vector<MLine*> myLines;// MLine's of the current nonloop GEdge
+ std::set<MVertex*>::iterator itf = first.begin();
+ MVertex* cv0 = *itf;// starting with a first vertex
+ while(last.find(cv0) != last.end()){// until reaching the corresponding last vertex
+ myLines.push_back(v02edg[cv0]);// push the correspong MLine
+ v02edg.erase(cv0);// and erasing it from the "list"
+ cv0 = myLines.back()->getVertex(1);// next vertex
+ }
+
+ v0 = this->model()->addVertex((*itf)->x(),myLines[0]->getVertex(0)->y(),myLines[0]->getVertex(0)->z(),1.);
+ v1 = this->model()->addVertex(cv0->x(),cv0->y(),cv0->z(),1.);
+
+ GEdge* toadd = this->model()->addLine(v0,v1);// new GEdge
+ toadd->lines = myLines;// associated MLine's
+
+ partition[i]->my_GEdges.push_back(toadd);
+ partition[ii]->my_GEdges.push_back(toadd);
+
+ first.erase(itf);// next first vertex of a nonloop GEdge
+ }//end while first.empty()
+
+ std::vector<MLine*> my_MLines;// remaining MLines for 'loop'GEdge's
+ MVertex* cv0 = v02edg.begin()->first;
+ while(v02edg.find(cv0) != v02edg.end()){// a loop GEdge
+ my_MLines.push_back(v02edg[cv0]);// current MLine of the loop
+ v02edg.erase(cv0);// update the container
+ cv0 = my_MLines.back()->getVertex(1);// next MLine of the loop
+ }
+
+ v0 = this->model()->addVertex(cv0->x(),cv0->y(),cv0->z(),1.);
+ v1 = v0;
+
+ GEdge* to_add = this->model()->addLine(v0,v1);
+ to_add->lines = my_MLines;
+
+ partition[i]->my_GEdges.push_back(to_add);
+ partition[ii]->my_GEdges.push_back(to_add);
+
+ }// end while v02edg.empty()
+
+ }//end for ii
+
+ }// end for i
return;
}
diff --git a/Geo/discreteFace.h b/Geo/discreteFace.h
index 8b4abfb2d473e9ee83585593c8760afdab9f8431..5b7071fdd33c34a12a7e19d86794fa4e0e848c35 100644
--- a/Geo/discreteFace.h
+++ b/Geo/discreteFace.h
@@ -15,6 +15,7 @@
#include "meshPartition.h"
#include "MTriangle.h"
#include "MEdge.h"
+#include "MLine.h"
#include <stdlib.h>
class discreteDiskFace;