From a558f142ae0e4e4b50bcc4719984be30abd0d6a0 Mon Sep 17 00:00:00 2001
From: Pierre-Alexandre Beaufort <pierre-alexandre.beaufort@uclouvain.be>
Date: Wed, 22 Jun 2016 11:41:22 +0000
Subject: [PATCH] Splitting for large aspect ratio
---
Geo/discreteDiskFace.cpp | 54 +++----------
Geo/discreteDiskFace.h | 34 ++++++--
Geo/discreteFace.cpp | 170 ++++++++++++++++++++++++++++++++-------
Geo/discreteFace.h | 4 +
4 files changed, 183 insertions(+), 79 deletions(-)
diff --git a/Geo/discreteDiskFace.cpp b/Geo/discreteDiskFace.cpp
index eaf94c0f2d..7094c8fa19 100644
--- a/Geo/discreteDiskFace.cpp
+++ b/Geo/discreteDiskFace.cpp
@@ -220,16 +220,12 @@ discreteDiskFace::discreteDiskFace(GFace *gf, triangulation* diskTriangulation,
_order = p;
_N = (p+1)*(p+2)/2;
discrete_triangles.resize(mesh.size());
-
std::map<MVertex*,MVertex*> v2v;// mesh vertex |-> face vertex
std::map<MEdge,MVertex*,Less_Edge> ed2nodes; // edge to interior node(s)
for (unsigned int i=0;i<mesh.size();i++){ // triangle by triangle
-
std::vector<MVertex*> vs; // MTriangle vertices
for (unsigned int j=0; j<3; j++){ // loop over vertices AND edges of the current triangle
-
MVertex *v = mesh[i]->getVertex(j);// firstly, edge vertices
-
if (v->onWhat()->dim() == 2) {
std::map<MVertex*,MVertex*> :: iterator it = v2v.find(v);
if (it == v2v.end()){
@@ -255,8 +251,7 @@ discreteDiskFace::discreteDiskFace(GFace *gf, triangulation* diskTriangulation,
}
else vs.push_back(v);
}
- else vs.push_back(v);
-
+ else vs.push_back(v);
}
if (_order == 2) {// then, interior nodes :-)
for (unsigned int ie=0; ie<3; ie++){ // firstly, edge interior nodes :-)
@@ -269,25 +264,20 @@ discreteDiskFace::discreteDiskFace(GFace *gf, triangulation* diskTriangulation,
discrete_vertices.push_back(mv);
ed2nodes[me]=mv;
}
- else{
- vs.push_back(it->second);
- }
+ else vs.push_back(it->second);
}
}// end order == 2
if (_order==1)
discrete_triangles[i] = new MTriangle (vs);
else if (_order==2)
discrete_triangles[i] = new MTriangle6 (vs);
+
}// end loop over triangles
-
geoTriangulation = new triangulation(discrete_triangles,gf);
-
allNodes = geoTriangulation->vert;
_totLength = geoTriangulation->bord.rbegin()->first;
_U0 = geoTriangulation->bord.rbegin()->second;
-
orderVertices(_totLength, _U0, _coords);
-
parametrize(false);
buildOct(CAD);
@@ -295,12 +285,12 @@ discreteDiskFace::discreteDiskFace(GFace *gf, triangulation* diskTriangulation,
Msg::Info("discreteDiskFace:: parametrization is not one-to-one; fixing "
"the discrete system.");
parametrize(true);
- buildOct(CAD);
- if(!checkOrientationUV()) Msg::Fatal("discreteDiskFace:: failing to fix the discrete system");
+ buildOct(CAD);
}
- putOnView(true,true);
+ putOnView(true,false);
printParamMesh();
+ if(!checkOrientationUV()) Msg::Fatal("discreteDiskFace:: failing to fix the discrete system");
}
/*end BUILDER*/
@@ -364,11 +354,11 @@ bool discreteDiskFace::parametrize(bool one2one) const
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){
@@ -413,7 +403,6 @@ bool discreteDiskFace::parametrize(bool one2one) const
lsys_u->systemSolve();
lsys_v->systemSolve();
Msg::Debug("Systems solved");
-
for(std::set<MVertex *>::iterator itv = allNodes.begin(); itv !=allNodes.end() ; ++itv){
MVertex *v = *itv;
double value_U, value_V;
@@ -443,10 +432,7 @@ void discreteDiskFace::getTriangleUV(const double u,const double v,
double &_xi, double &_eta)const{
double uv[3] = {u,v,0.};
*mt = (discreteDiskFaceTriangle*) Octree_Search(uv,oct);
- if (!(*mt)){
- //Msg::Error("discreteDiskFace::getTriangleUV(), the octree does not have a triangle containing (u;v)=(%f;%f)",u,v);
- return;
- }
+ if (!(*mt)) return;
double Xi[2];
double U[2] = {u,v};
@@ -478,6 +464,7 @@ bool discreteDiskFace::checkOrientationUV()
p3[0] = ct->p[2].x(); p3[1] = ct->p[2].y();
current = robustPredicates::orient2d(p1, p2, p3);
if(initial*current < 0.) {
+ Msg::Error("Triangle UV %d has not the correct orientation",i+1);
return false;
break;
}
@@ -787,27 +774,6 @@ void discreteDiskFace::putOnView(bool Xu, bool Ux)
fclose(UVz);
}
}
- /*
- #ifdef HAVE_POST
- std::map<int, std::vector<double> > u;
- std::map<int, std::vector<double> > v;
- for(std::set<MVertex*>::iterator iv = allNodes.begin(); iv!=allNodes.end(); ++iv){
- MVertex *p = *iv;
- u[p->getNum()].push_back(coordinates[p].x());
- v[p->getNum()].push_back(coordinates[p].y());
- }
- PView* view_u = new PView("u", "NodeData", GModel::current(), u);
- PView* view_v = new PView("v", "NodeData", GModel::current(), v);
- view_u->setChanged(true);
- view_v->setChanged(true);
- view_u->write("param_u.msh", 5);
- view_v->write("param_v.msh", 5);
- delete view_u;
- delete view_v;
- #else
- Msg::Error("discreteDiskFace: cannot export without post module")
- #endif
- */
}
// useful for mesh generators
diff --git a/Geo/discreteDiskFace.h b/Geo/discreteDiskFace.h
index 23a01ebf4d..75292a1be1 100644
--- a/Geo/discreteDiskFace.h
+++ b/Geo/discreteDiskFace.h
@@ -48,6 +48,8 @@ class triangulation {
public:
// attributes
+ double area;
+ bool seamPoint;
std::vector<MElement*> tri;// triangles
std::set<MVertex*> vert;// nodes
// edge to 1 or 2 triangle(s), their num into the vector of MElement*
@@ -59,7 +61,14 @@ class triangulation {
std::list<GEdge*> my_GEdges;
- //methods
+ //---- methods
+
+ double aspectRatio()
+ {
+ double L = bord.rbegin()->first;
+ return L*L / (area*M_PI);
+ }
+
int genus()
{
return ( ed2tri.size() - vert.size() - tri.size() + 2 - bord.size() )/2;
@@ -69,11 +78,16 @@ class triangulation {
{
for(unsigned int i = 0; i < tri.size(); ++i){
MElement* t = tri[i];
+ SPoint3 P[3];
for(int j = 0; j < t->getNumVertices() ; j++){
MVertex* tv = t->getVertex(j);
+ if(j<3)
+ P[j].setPosition(tv->x(),tv->y(),tv->z());
std::set<MVertex*>::iterator it = vert.find(tv);
if (it == vert.end()) vert.insert (tv);
}
+ SVector3 N(crossprod(P[1]-P[0],P[2]-P[0]));
+ area += N.norm();
}
}
@@ -114,13 +128,16 @@ class triangulation {
vecver.erase(vecver.begin());
std::map<MVertex*,std::vector<MVertex*> >::iterator im = firstNode2Edge.find(first);
- if (im != firstNode2Edge.end())
- Msg::Error("Incorrect topology in discreteDiskFace %d", gf->tag());
+ if (im != firstNode2Edge.end()){
+ Msg::Info("Supposed seam point for discreteFace %d", gf->tag());
+ seamPoint = true;
+ break;
+ }
firstNode2Edge[first] = vecver;
firstNode2Edge[first].push_back(last);
}
- while (!firstNode2Edge.empty()) {
+ while (!firstNode2Edge.empty()&&!seamPoint) {
std::vector<MVertex*> loop;
double length = 0.;
@@ -148,6 +165,8 @@ class triangulation {
void assign()
{
+ area = 0.;
+ seamPoint = false;
assignVert();
assignEd2tri();
assignBord();
@@ -194,11 +213,14 @@ class discreteDiskFace : public GFace {
GPoint intersectionWithCircle(const SVector3 &n1, const SVector3 &n2,
const SVector3 &p, const double &d,
double uv[2]) const;
+
+
+ std::vector<MElement*> discrete_triangles;
protected:
- // a copy of the mesh that should not be destroyed
+ // a copy of the mesh that should not be destroyed
triangulation* initialTriangulation;
triangulation* geoTriangulation;// parametrized triangulation
- std::vector<MElement*> discrete_triangles;
+
std::vector<MVertex*> discrete_vertices;
int _order;
diff --git a/Geo/discreteFace.cpp b/Geo/discreteFace.cpp
index e42a9414c4..e270b06d75 100644
--- a/Geo/discreteFace.cpp
+++ b/Geo/discreteFace.cpp
@@ -20,6 +20,8 @@ extern "C" {
}
#endif
+int cp=0;
+
discreteFace::discreteFace(GModel *model, int num) : GFace(model, num)
{
Surface *s = Create_Surface(num, MSH_SURF_DISCRETE);
@@ -110,9 +112,9 @@ void discreteFace::secondDer(const SPoint2 ¶m,
void discreteFace::createGeometry()
{
checkAndFixOrientation();
-
int order = 1;
int nPart = 2;
+ double eta = .125;
#if defined(HAVE_SOLVER) && defined(HAVE_ANN)
@@ -124,13 +126,8 @@ void discreteFace::createGeometry()
std::vector<MElement*> tem(triangles.begin(),triangles.end());
triangulation* init = new triangulation(tem,this);
-
toSplit.push(init);
- //int mygen, compteur=1;//#debug
- //Msg::Info("First Genus Value:");
- //mygen=1; //#debug
- //if(mygen!=0){// #debug
- if((toSplit.top())->genus()!=0){
+ if((toSplit.top())->genus()!=0 || (toSplit.top())->aspectRatio() < eta || (toSplit.top())->seamPoint){
while( !toSplit.empty()){
std::vector<triangulation*> part;
@@ -138,19 +135,15 @@ void discreteFace::createGeometry()
toSplit.pop();
split(tosplit,part,nPart);
- delete tosplit; // #mark
+ delete tosplit;
for(unsigned int i=0; i<part.size(); i++){
- //Msg::Info("Partition %d Genus:",compteur);//#debug
- //std::cin>>mygen;//#debug
- //if (mygen!=0)//#debug
- if(part[i]->genus()!=0)
+ if(part[i]->genus()!=0 || part[i]->aspectRatio() < eta || part[i]->seamPoint)
toSplit.push(part[i]);
else{
toParam.push_back(part[i]);
part[i]->idNum=id++;
}
- //compteur++; //#debug
}// end for i
}// !.empty()
}// end if it is not disk-like
@@ -158,12 +151,11 @@ void discreteFace::createGeometry()
toParam.push_back(toSplit.top());
toSplit.top()->idNum=id++;
}
-
- updateTopology(toParam);
-
+ updateTopology(toParam);
for(unsigned int i=0; i<toParam.size(); i++){
+ fillHoles(toParam[i]);
std::vector<MElement*> mytri = toParam[i]->tri;
- discreteDiskFace *df = new discreteDiskFace (this,toParam[i], order,(_CAD.empty() ? NULL : &_CAD));
+ discreteDiskFace *df = new discreteDiskFace (this,toParam[i], order,(_CAD.empty() ? NULL : &_CAD));
df->replaceEdges(toParam[i]->my_GEdges);
_atlas.push_back(df);
}
@@ -181,7 +173,7 @@ void discreteFace::gatherMeshes()
for (unsigned int j=0;j<_atlas[i]->triangles.size(); j++){
MTriangle *t = _atlas[i]->triangles[j];
SPoint2 p0,p1,p2;
- reparamMeshVertexOnFace(t->getVertex(0),_atlas[i], p0);
+ reparamMeshVertexOnFace(t->getVertex(0),_atlas[i], p0);
reparamMeshVertexOnFace(t->getVertex(1),_atlas[i], p1);
reparamMeshVertexOnFace(t->getVertex(2),_atlas[i], p2);
SPoint2 pc = (p0+p1+p2)*(1./3.0);
@@ -189,8 +181,7 @@ void discreteFace::gatherMeshes()
double xi, eta;
_atlas[i]->getTriangleUV(pc.x(),pc.y(), &mt, xi,eta);
if (mt && mt->gf)mt->gf->triangles.push_back(t);
- else Msg::Warning ("FILE %s LINE %d Triangle has no classification",__FILE__,__LINE__);
-
+ else Msg::Warning ("FILE %s LINE %d Triangle from atlas part %u has no classification for (u;v)=(%f;%f)",__FILE__,__LINE__,i+1,pc.x(),pc.y());
}
@@ -218,14 +209,96 @@ void discreteFace::mesh(bool verbose)
{
#if defined(HAVE_ANN) && defined(HAVE_SOLVER) && defined(HAVE_MESH)
if (!CTX::instance()->meshDiscrete) return;
- for (unsigned int i=0;i<_atlas.size();i++)
+ for (unsigned int i=0;i<_atlas.size();i++){
+ Msg::Info("Discrete Face %d is going to be meshed",i);
+ printAtlasMesh(_atlas[i]->discrete_triangles,i);
_atlas[i]->mesh(verbose);
-
+ //printAtlasMesh(_atlas[i],i);
+ }
+
gatherMeshes();
meshStatistics.status = GFace::DONE;
#endif
}
+void discreteFace::printAtlasMesh(std::vector<MElement*> elm, int I)
+{
+
+ std::map<MVertex*,int> mv2int;
+ char buffer[16];
+ sprintf(buffer,"atlas_mesh%d.msh",I);
+ FILE* pmesh = Fopen(buffer,"w");
+
+ std::set<MVertex*> meshvertices;
+
+ for(unsigned int i=0; i<elm.size(); ++i){
+ MElement* tri = elm[i];
+ for(unsigned int j=0; j<3; j++)
+ if (meshvertices.find(tri->getVertex(j))==meshvertices.end()) meshvertices.insert(tri->getVertex(j));
+ }
+
+ fprintf(pmesh,"$MeshFormat\n2.2 0 8\n$EndMeshFormat\n$Nodes\n%u\n",(unsigned int)meshvertices.size());
+ int count = 1;
+ for(std::set<MVertex*>::iterator it = meshvertices.begin(); it!=meshvertices.end(); ++it){
+ fprintf(pmesh,"%d %f %f %f\n",count,(*it)->x(),(*it)->y(),(*it)->z());
+ mv2int[*it] = count;
+ count++;
+ }
+ fprintf(pmesh,"$EndNodes\n$Elements\n%u\n",elm.size());
+ for(unsigned int i=0; i<elm.size(); i++){
+ MElement* tri = elm[i];
+ fprintf(pmesh,"%d 2 2 0 0",i+1);
+ for(int j=0; j<3; j++){
+ MVertex* mv = tri->getVertex(j);
+ fprintf(pmesh," %d",mv2int[mv]);
+ }
+ fprintf(pmesh,"\n");
+ }
+ fprintf(pmesh,"$EndElements\n");
+ fclose(pmesh);
+
+}
+
+
+
+void discreteFace::printAtlasMesh(discreteDiskFace* ddf, int I)
+{
+
+ std::map<MVertex*,int> mv2int;
+ char buffer[16];
+ sprintf(buffer,"atlas_mesh%d.msh",I);
+ FILE* pmesh = Fopen(buffer,"w");
+
+ std::set<MVertex*> meshvertices;
+
+ for(unsigned int i=0; i<ddf->triangles.size(); ++i){
+ MTriangle* tri = ddf->triangles[i];
+ for(unsigned int j=0; j<3; j++)
+ if (meshvertices.find(tri->getVertex(j))==meshvertices.end()) meshvertices.insert(tri->getVertex(j));
+ }
+
+ fprintf(pmesh,"$MeshFormat\n2.2 0 8\n$EndMeshFormat\n$Nodes\n%u\n",(unsigned int)meshvertices.size());
+ int count = 1;
+ for(std::set<MVertex*>::iterator it = meshvertices.begin(); it!=meshvertices.end(); ++it){
+ fprintf(pmesh,"%d %f %f %f\n",count,(*it)->x(),(*it)->y(),(*it)->z());
+ mv2int[*it] = count;
+ count++;
+ }
+ fprintf(pmesh,"$EndNodes\n$Elements\n%u\n",ddf->triangles.size());
+ for(unsigned int i=0; i<ddf->triangles.size(); i++){
+ MTriangle* tri = ddf->triangles[i];
+ fprintf(pmesh,"%d 2 2 0 0",i+1);
+ for(int j=0; j<3; j++){
+ MVertex* mv = tri->getVertex(j);
+ fprintf(pmesh," %d",mv2int[mv]);
+ }
+ fprintf(pmesh,"\n");
+ }
+ fprintf(pmesh,"$EndElements\n");
+ fclose(pmesh);
+
+}
+
void discreteFace::checkAndFixOrientation(){
@@ -296,7 +369,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("discreteFace: triangle %d has been reoriented.",neigs[i]->getNum());
+ //Msg::Info("discreteFace: triangle %d has been reoriented.",neigs[i]->getNum());
}
break;
}
@@ -323,7 +396,7 @@ void discreteFace::setupDiscreteEdge(discreteEdge*de,std::vector<MLine*>mlines,s
de->mesh_vertices.push_back(mlines[i]->getVertex(0));
if(trash) trash->insert(mlines[i]->getVertex(0));
}
- de->createGeometry();
+ de->createGeometry();
}
@@ -380,7 +453,10 @@ void discreteFace::splitDiscreteEdge ( GEdge *de , GVertex *gv, discreteEdge* ne
void discreteFace::split(triangulation* trian,std::vector<triangulation*> &partition,int nPartitions)
{
#if defined(HAVE_SOLVER) && defined(HAVE_ANN) && defined(HAVE_METIS)
-
+ /*
+ if(trian->aspectRatio() < 4.)
+ printf("it is splitted because of aspect ratio");
+ */
int nVertex = trian->tri.size(); // number of elements
int nEdge = trian->ed2tri.size() - trian->borderEdg.size();// number of edges, (without the boundary ones)
@@ -422,8 +498,9 @@ void discreteFace::split(triangulation* trian,std::vector<triangulation*> &parti
elem[part[i]].push_back(trian->tri[i]);
el2part[trian->tri[i]] = part[i];
}
+
//check connectivity
- for(int p=0; p<nPartitions; p++){// part by part
+ for(int p=0; p<nPartitions; p++){// part by part
std::set<MElement*> celem(elem[p].begin(),elem[p].end());// current elements of the p-th part
std::queue<MElement*> my_todo; // todo list, for adjacency check - in order to check the connectivity of the part
std::map<MElement*,bool> check_todo; // help to complete todo list
@@ -449,7 +526,7 @@ void discreteFace::split(triangulation* trian,std::vector<triangulation*> &parti
}// end for j
}// end while
if(!celem.empty()){// if the set is empty, it means that the part was connected
- Msg::Info("discreteFace::split(), a partition is not connected: it is going to be fixed.");
+ Msg::Info("discreteFace::split(), a partition (%d) is not connected: it is going to be fixed.",p);
std::vector<MElement*> relem(celem.begin(),celem.end());// new part
for(unsigned int ie=0; ie<relem.size(); ie++)// updating the id part of the element belonging to the new part
el2part[relem[ie]] = nPartitions;
@@ -463,8 +540,6 @@ void discreteFace::split(triangulation* trian,std::vector<triangulation*> &parti
elem[p] = upe;
}
}// end for p
-
-
for(int i=0; i<nPartitions; i++)// new triangulation of the connected parts
partition.push_back(new triangulation(elem[i],this));
@@ -603,6 +678,43 @@ void discreteFace::updateTopology(std::vector<triangulation*>&partition)
#endif
}
+
+void discreteFace::fillHoles(triangulation* trian)
+{
+
+ std::map<double,std::vector<MVertex*> > bords = trian->bord;
+ std::map<double,std::vector<MVertex*> >::reverse_iterator it = bords.rbegin();
+ ++it;
+ for(; it!=bords.rend(); ++it){
+ double x[3] = {0.,0.,0.};
+ std::vector<MVertex*> mv = it->second;
+ for(unsigned int j=0; j<mv.size(); j++){
+ x[0] += mv[j]->x();
+ x[1] += mv[j]->y();
+ x[2] += mv[j]->z();
+ }
+ x[0] /= mv.size(); x[1] /= mv.size(); x[2] /= mv.size();
+ MVertex* center = new MVertex(x[0],x[1],x[2]);
+ this->mesh_vertices.push_back(center);
+ center->setEntity(this);
+ trian->vert.insert(center);
+ for(unsigned int j=1; j<mv.size(); j++)
+ addTriangle(trian,new MTriangle(mv[j],mv[j-1],center));
+ addTriangle(trian,new MTriangle(mv[0],mv[mv.size()-1],center));
+ }
+
+}
+
+void discreteFace::addTriangle(triangulation* trian, MTriangle* t)
+{// #mark quid borders ?
+ for(int i=0; i<3; i++){
+ MEdge ed = t->getEdge(i);
+ trian->ed2tri[ed].push_back(trian->tri.size());
+ }
+ trian->tri.push_back(t);
+}
+
+
// delete all discrete disk faces
//void discreteFace::deleteAtlas() {
//}
diff --git a/Geo/discreteFace.h b/Geo/discreteFace.h
index d043ffc2ba..1f6b0f61af 100644
--- a/Geo/discreteFace.h
+++ b/Geo/discreteFace.h
@@ -36,6 +36,8 @@ class discreteFace : public GFace {
void splitDiscreteEdge(GEdge*,GVertex*,discreteEdge*[2]);
void updateTopology(std::vector<triangulation*>&);
void split(triangulation*,std::vector<triangulation*>&,int);
+ void fillHoles(triangulation*);
+ void addTriangle(triangulation*,MTriangle*);
GPoint point(double par1, double par2) const;
SPoint2 parFromPoint(const SPoint3 &p, bool onSurface=true) const;
SVector3 normal(const SPoint2 ¶m) const;
@@ -53,6 +55,8 @@ class discreteFace : public GFace {
void createGeometry();
void gatherMeshes();
virtual void mesh (bool verbose);
+ void printAtlasMesh(std::vector<MElement*>, int);
+ void printAtlasMesh(discreteDiskFace*, int);
std::vector<discreteDiskFace*> _atlas;
std::vector<GFace*> _CAD;
};
--
GitLab