From d95bf7f8f1ff2b5f5d584dfc40f65fa835235d0c Mon Sep 17 00:00:00 2001
From: Pierre-Alexandre Beaufort <pierre-alexandre.beaufort@uclouvain.be>
Date: Tue, 19 Jul 2016 13:51:35 +0000
Subject: [PATCH] (likely) robust filler of holes
---
Geo/discreteDiskFace.cpp | 37 ++---
Geo/discreteDiskFace.h | 34 ++++-
Geo/discreteFace.cpp | 300 ++++++++++++++++++++++++++++++++++++++-
Geo/discreteFace.h | 7 +
4 files changed, 352 insertions(+), 26 deletions(-)
diff --git a/Geo/discreteDiskFace.cpp b/Geo/discreteDiskFace.cpp
index 8bc019a514..4c28872ce9 100644
--- a/Geo/discreteDiskFace.cpp
+++ b/Geo/discreteDiskFace.cpp
@@ -274,6 +274,7 @@ discreteDiskFace::discreteDiskFace(int id, GFace *gf, triangulation* diskTriangu
}// end loop over triangles
geoTriangulation = new triangulation(id,discrete_triangles,gf);
+ geoTriangulation->fillingHoles = diskTriangulation->fillingHoles;
allNodes = geoTriangulation->vert;
_totLength = geoTriangulation->bord.rbegin()->first;
_U0 = geoTriangulation->bord.rbegin()->second;
@@ -281,6 +282,7 @@ discreteDiskFace::discreteDiskFace(int id, GFace *gf, triangulation* diskTriangu
parametrize(false);
buildOct(CAD);
+
if (!checkOrientationUV()){
Msg::Info("discreteDiskFace:: parametrization is not one-to-one; fixing "
"the discrete system.");
@@ -312,21 +314,21 @@ void discreteDiskFace::buildOct(std::vector<GFace*> *CAD) const
oct = Octree_Create(maxElePerBucket, origin, ssize, discreteDiskFaceBB,
discreteDiskFaceCentroid, discreteDiskFaceInEle);
- _ddft = new discreteDiskFaceTriangle[discrete_triangles.size()];
-
- // if (CAD) printf("-------------> %ld %ld\n",CAD->size(),discrete_triangles.size());
-
+ _ddft = new discreteDiskFaceTriangle[discrete_triangles.size()-geoTriangulation->fillingHoles.size()];
+ int c = 0;
for(unsigned int i = 0; i < discrete_triangles.size(); ++i){
- MElement *t = discrete_triangles[i];
- discreteDiskFaceTriangle* my_ddft = &_ddft[i];
- my_ddft->p.resize(_N);
- for(int io=0; io<_N; io++)
- my_ddft->p[io] = coordinates[t->getVertex(io)];
+ if(geoTriangulation->fillingHoles.find(i)==geoTriangulation->fillingHoles.end()){
+ MElement *t = discrete_triangles[i];
+ discreteDiskFaceTriangle* my_ddft = &_ddft[c++];
+ my_ddft->p.resize(_N);
+ for(int io=0; io<_N; io++)
+ my_ddft->p[io] = coordinates[t->getVertex(io)];
- my_ddft->gf = CAD ? (*CAD)[i] : _parent;
- my_ddft->tri = t;
+ my_ddft->gf = CAD ? (*CAD)[i] : _parent;
+ my_ddft->tri = t;
- Octree_Insert(my_ddft, oct);
+ Octree_Insert(my_ddft, oct);
+ }
}
Octree_Arrange(oct);
}
@@ -460,7 +462,7 @@ bool discreteDiskFace::checkOrientationUV()
double p3[2] = {ct->p[2].x(), ct->p[2].y()};
initial = robustPredicates::orient2d(p1, p2, p3);
unsigned int i=1;
- for (; i<discrete_triangles.size(); i++){
+ for (; i<discrete_triangles.size()-geoTriangulation->fillingHoles.size(); i++){
ct = &_ddft[i];
p1[0] = ct->p[0].x(); p1[1] = ct->p[0].y();
p2[0] = ct->p[1].x(); p2[1] = ct->p[1].y();
@@ -478,7 +480,7 @@ bool discreteDiskFace::checkOrientationUV()
double min, max;
std::vector<MVertex*> localVertices;
localVertices.resize(_N);
- for(unsigned int i=0; i<discrete_triangles.size(); i++){
+ for(unsigned int i=0; i<discrete_triangles.size()-geoTriangulation->fillingHoles.size(); i++){
ct = &_ddft[i];
for(int j=0; j<_N; j++)
localVertices[j] = new MVertex(ct->p[j].x(),ct->p[j].y(),0.);
@@ -689,7 +691,7 @@ void discreteDiskFace::putOnView(bool Xu, bool Ux)
fprintf(UVy,"View \"y(U)\"{\n");
fprintf(UVz,"View \"z(U)\"{\n");
}
- for (unsigned int i=0; i<discrete_triangles.size(); i++){
+ for (unsigned int i=0; i<discrete_triangles.size()-geoTriangulation->fillingHoles.size(); i++){
discreteDiskFaceTriangle* my_ddft = &_ddft[i];
if (_order == 1){
if(Xu){
@@ -787,7 +789,7 @@ GPoint discreteDiskFace::intersectionWithCircle(const SVector3 &n1, const SVecto
{
SVector3 n = crossprod(n1,n2);
n.normalize();
- for (int i=-1;i<(int)discrete_triangles.size();i++){
+ for (int i=-1;i<(int)(discrete_triangles.size()-geoTriangulation->fillingHoles.size());i++){
discreteDiskFaceTriangle *ct = NULL;
double U,V;
if (i == -1) getTriangleUV(uv[0],uv[1], &ct, U,V);
@@ -938,6 +940,7 @@ static MEdge getEdge (MElement *t, MVertex *v){
if (t->getVertex(i) == v) return t->getEdge((i+1)%3);
}
+/* warning
static double computeDistanceLinePoint (MVertex *v1, MVertex *v2, MVertex *v){
SVector3 U = v2->point() - v1->point();
@@ -947,7 +950,7 @@ static double computeDistanceLinePoint (MVertex *v1, MVertex *v2, MVertex *v){
return xx.norm() / U.norm();
}
-
+*/
static double computeDistance (MVertex *v1, double d1, MVertex *v2, double d2, MVertex *v){
diff --git a/Geo/discreteDiskFace.h b/Geo/discreteDiskFace.h
index c8fc29dc31..d5e3949045 100644
--- a/Geo/discreteDiskFace.h
+++ b/Geo/discreteDiskFace.h
@@ -23,6 +23,20 @@
#include "PView.h"
#include "robustPredicates.h"
+inline double tri3Darea(MVertex* mv0, MVertex* mv1, MVertex* mv2)
+{
+
+ SVector3 v1(mv1->x()-mv0->x(),mv1->y()-mv0->y(),mv1->z()-mv0->z());
+ SVector3 v2(mv2->x()-mv0->x(),mv2->y()-mv0->y(),mv2->z()-mv0->z());
+
+ SVector3 n(v1.y()*v2.z()-v2.y()*v1.z(),v2.x()*v1.z()-v1.x()*v2.z(),v1.x()*v2.y()-v2.x()*v1.y());
+
+
+ return .5*n.norm();
+
+}
+
+
inline int nodeLocalNum(MElement* e, MVertex* v)
{
for(int i=0; i<e->getNumVertices(); i++)
@@ -39,6 +53,15 @@ inline int edgeLocalNum(MElement* e, MEdge ed)
return -1;
}
+inline MEdge maxEdge(MElement* e)
+{
+ MEdge maxEd = e->getEdge(0);
+ for(int i=0; i<e->getNumEdges(); i++)
+ if (maxEd.length() < e->getEdge(i).length())
+ maxEd = e->getEdge(i);
+ return maxEd;
+}
+
class ANNkd_tree;
class Octree;
class GRbf;
@@ -48,22 +71,23 @@ class triangulation {
public:
// attributes
+ int idNum; // number of identification, for hashing purposes
+ std::vector<MElement*> tri;// trianglse
+ GFace* gf;
+
double maxD;
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*
std::map<MEdge,std::vector<int>,Less_Edge> ed2tri;
std::map<double,std::vector<MVertex*> > bord; //border(s)
std::set<MEdge,Less_Edge> borderEdg; // border edges
- GFace *gf;
- int idNum; // number of identification, for hashing purposes
-
+
std::list<GEdge*> my_GEdges;
+ std::set<int> fillingHoles;
//---- methods
-
double geodesicDistance ();
double aspectRatio()
diff --git a/Geo/discreteFace.cpp b/Geo/discreteFace.cpp
index 22e70ab8e1..a6e6196cec 100644
--- a/Geo/discreteFace.cpp
+++ b/Geo/discreteFace.cpp
@@ -20,7 +20,6 @@ extern "C" {
}
#endif
-int cp=0;
discreteFace::discreteFace(GModel *model, int num) : GFace(model, num)
{
@@ -216,6 +215,7 @@ void discreteFace::gatherMeshes()
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++){
@@ -707,20 +707,312 @@ void discreteFace::fillHoles(triangulation* trian)
this->mesh_vertices.push_back(center);
center->setEntity(this);
trian->vert.insert(center);
- for(unsigned int j=1; j<mv.size(); j++)
+ for(unsigned int j=1; j<mv.size(); j++){
+ if(std::abs(tri3Darea(mv[j],mv[j-1],center))<1e-15){
+ MVertex temp(center->x()+mv[j]->x()/mv.size(),center->y()+mv[j-1]->y()/mv.size(),(1+1./mv.size())*center->z());
+ *center = temp;
+ center->setEntity(this);
+ }
addTriangle(trian,new MTriangle(mv[j],mv[j-1],center));
+ }
+ if(std::abs(tri3Darea(mv[0],mv[mv.size()-1],center))<1e-15){
+ MVertex temp(center->x()+mv[0]->x()/mv.size(),center->y()+mv[mv.size()-1]->y()/mv.size(),(1+1./mv.size())*center->z());
+ *center = temp;
+ center->setEntity(this);
+ }
addTriangle(trian,new MTriangle(mv[0],mv[mv.size()-1],center));
}
#endif
}
+void discreteFace::bisectionEdg(triangulation* trian)
+{
+ std::map<MEdge,double,Less_Edge> ed2angle;
+ std::map<double,std::vector<MEdge> > ed2refine;
+ /*
+ checkEdges
+ for each ed in ed2refine
+ ...while ed exists
+ ......terminal from lepp(ed)
+ ......refine terminal
+ ......update ed2refine
+ */
+ int count = 1;
+ checkEdgesBis(trian,ed2refine,ed2angle);
+ while(!ed2refine.empty()){
+ //for(it = ed2refine.rbegin(); it != ed2refine.rend(); ++it){
+ std::map<double,std::vector<MEdge> >::reverse_iterator it = ed2refine.rbegin();
+ double key = it->first;
+ std::vector<MEdge> & ved = it->second;
+ while(!ved.empty()){
+ //for(unsigned int i=0; i<ved.size(); i++){
+ MEdge ed = ved[0];
+ while(ed2angle.find(ed) != ed2angle.end()){
+ MEdge terminal = lepp(trian,ed);
+ ed2angle.erase(terminal);
+ std::vector<MEdge> ed2update;
+ refineTriangles(trian,terminal,ed2update);
+ updateEd2refineBis(trian,ed2angle,ed2update,ed2refine);
+ printAtlasMesh(trian->tri,-count++);
+ }// end while exists
+ ed2refine[key].erase(ved.begin());
+ ved = ed2refine[key];
+ }// end for unsigned
+ ed2refine.erase(key);
+ }// end for it
+
+}
+
+bool discreteFace::checkEdges(triangulation* trian,std::map<double,std::vector<MEdge> >& ed2refine,std::map<MEdge,double,Less_Edge>& ed2angle)
+{
+
+ for(unsigned int i=0; i<trian->tri.size(); i++){
+ MElement* t = trian->tri[i];
+ std::vector<MEdge> eds;
+ eds.push_back(t->getEdge(0));eds.push_back(t->getEdge(1));eds.push_back(t->getEdge(2));
+ for(int j=0; j<3; j++){
+ if(trian->ed2tri[eds[j]].size()>1){
+ MEdge ed1 = eds[j+1>2?j-2:j+1];
+ MEdge ed2 = eds[j+2>2?j-1:j+2];
+ SVector3 v1(ed1.getVertex(0)->x()-ed1.getVertex(1)->x(),
+ ed1.getVertex(0)->y()-ed1.getVertex(1)->y(),
+ ed1.getVertex(0)->z()-ed1.getVertex(1)->z());
+ SVector3 v2(ed2.getVertex(1)->x()-ed2.getVertex(0)->x(),
+ ed2.getVertex(1)->y()-ed2.getVertex(0)->y(),
+ ed2.getVertex(1)->z()-ed2.getVertex(0)->z());
+ double dotp = dot(v1,v2);
+ double lp = ed1.length()*ed2.length();
+ std::map<MEdge,double,Less_Edge>::iterator check = ed2angle.find(eds[j]);
+ if(check!=ed2angle.end()){
+ ed2angle[eds[j]] += std::acos(dotp/lp);
+ if(ed2angle[eds[j]]>M_PI)
+ ed2refine[eds[j].length()].push_back(eds[j]);
+ else
+ ed2angle.erase(eds[j]);
+ }
+ else
+ ed2angle[eds[j]] = std::acos(dotp/lp);
+ }//end if size()>1
+ }// end for j
+ }//end for tri
+ return !(ed2refine.empty());
+}
+
+
+bool discreteFace::checkEdgesBis(triangulation* trian,std::map<double,std::vector<MEdge> >& ed2refine,std::map<MEdge,double,Less_Edge>& ed2angle)
+{
+
+ for(unsigned int i=0; i<trian->tri.size(); i++){
+ MElement* t = trian->tri[i];
+ std::vector<MEdge> eds;
+ eds.push_back(t->getEdge(0));eds.push_back(t->getEdge(1));eds.push_back(t->getEdge(2));
+ for(int j=0; j<3; j++){
+ MEdge ed1 = eds[j+1>2?j-2:j+1];
+ MEdge ed2 = eds[j+2>2?j-1:j+2];
+ SVector3 v1(ed1.getVertex(0)->x()-ed1.getVertex(1)->x(),
+ ed1.getVertex(0)->y()-ed1.getVertex(1)->y(),
+ ed1.getVertex(0)->z()-ed1.getVertex(1)->z());
+ SVector3 v2(ed2.getVertex(1)->x()-ed2.getVertex(0)->x(),
+ ed2.getVertex(1)->y()-ed2.getVertex(0)->y(),
+ ed2.getVertex(1)->z()-ed2.getVertex(0)->z());
+ double dotp = dot(v1,v2);
+ double lp = ed1.length()*ed2.length();
+ std::map<MEdge,double,Less_Edge>::iterator check = ed2angle.find(eds[j]);
+ if(check==ed2angle.end()){
+ ed2angle[eds[j]] = std::acos(dotp/lp);
+ if(ed2angle[eds[j]]>M_PI/2.)
+ ed2refine[eds[j].length()].push_back(eds[j]);
+ else
+ ed2angle.erase(eds[j]);
+ }
+ }// end for j
+ }//end for tri
+ return !(ed2refine.empty());
+}
+
+
+MEdge discreteFace::lepp(triangulation* trian,MEdge ed){
+ std::vector<int> intri = trian->ed2tri[ed];
+ bool isTerminal = false;
+ while(intri.size()>1 && !isTerminal){
+ MEdge ed0 = maxEdge(trian->tri[intri[0]]);
+ MEdge ed1 = maxEdge(trian->tri[intri[1]]);
+ if(ed != ed0){
+ ed = ed0;
+ intri = trian->ed2tri[ed];
+ }
+ else if(ed != ed1){
+ ed = ed1;
+ intri = trian->ed2tri[ed];
+ }
+ else isTerminal = true;
+ }// end while
+ return ed;
+}
+
+void discreteFace::refineTriangles(triangulation* trian,MEdge ed,std::vector<MEdge>&ed2update){
+
+ SPoint3 mid = ed.barycenter();
+ MVertex* mv = new MVertex(mid.x(),mid.y(),mid.z());
+ trian->vert.insert(mv);
+
+ std::vector<int> t = trian->ed2tri[ed];
+
+ MTriangle *t00, *t01;
+ int ie;
+
+ ie = edgeLocalNum(trian->tri[t[0]],ed);
+ t00 = new MTriangle(trian->tri[t[0]]->getEdge(ie+1 > 2 ? ie-2 : ie+1).getVertex(0),trian->tri[t[0]]->getEdge(ie+1 > 2 ? ie-2 : ie+1).getVertex(1),mv);
+ t01 = new MTriangle(trian->tri[t[0]]->getEdge(ie+2 > 2 ? ie-1 : ie+2).getVertex(0),trian->tri[t[0]]->getEdge(ie+2 > 2 ? ie-1 : ie+2).getVertex(1),mv);
+
+ trian->tri[t[0]] = t00;
+ trian->ed2tri[t00->getEdge(2)] = trian->ed2tri[ed];
+ trian->ed2tri[t00->getEdge(1)].push_back(t[0]);
+ trian->ed2tri[t00->getEdge(1)].push_back(trian->tri.size());
+
+ std::vector<int> oldtri = trian->ed2tri[ed];
+ oldtri[0] == t[0] ? oldtri[0] = trian->tri.size() : oldtri[1] = trian->tri.size();
+ trian->ed2tri[t01->getEdge(1)] = oldtri;
+ oldtri = trian->ed2tri[t01->getEdge(0)];
+ oldtri[0] == t[0] ? oldtri[0] = trian->tri.size() : oldtri[1] = trian->tri.size();
+ trian->ed2tri[t01->getEdge(0)] = oldtri;
+ trian->tri.push_back(t01);
+
+ ed2update.push_back(t00->getEdge(0));
+ ed2update.push_back(t01->getEdge(0));
+
+ if(t.size()>1){
+
+ ed2update.push_back(t01->getEdge(1));
+ ed2update.push_back(t00->getEdge(2));
+ MTriangle *t10, *t11;
+
+ this->mesh_vertices.push_back(mv);
+ mv->setEntity(this);
+
+ ie = edgeLocalNum(trian->tri[t[1]],ed);
+ t11 = new MTriangle(trian->tri[t[1]]->getEdge(ie+1 > 2 ? ie-2 : ie+1).getVertex(0),trian->tri[t[1]]->getEdge(ie+1 > 2 ? ie-2 : ie+1).getVertex(1),mv);
+ t10 = new MTriangle(trian->tri[t[1]]->getEdge(ie+2 > 2 ? ie-1 : ie+2).getVertex(0),trian->tri[t[1]]->getEdge(ie+2 > 2 ? ie-1 : ie+2).getVertex(1),mv);
+
+ trian->tri[t[1]] = t10;
+ trian->ed2tri[t10->getEdge(1)] = trian->ed2tri[ed];
+ trian->ed2tri[t10->getEdge(2)].push_back(t[1]);
+ trian->ed2tri[t10->getEdge(2)].push_back(trian->tri.size());
+
+ oldtri = trian->ed2tri[t11->getEdge(2)];
+ oldtri[0] == t[1] ? oldtri[0] = trian->tri.size() : oldtri[1] = trian->tri.size();
+ trian->ed2tri[t11->getEdge(2)] = oldtri;
+ oldtri = trian->ed2tri[t11->getEdge(0)];
+ oldtri[0] == t[1] ? oldtri[0] = trian->tri.size() : oldtri[1] = trian->tri.size();
+ trian->ed2tri[t11->getEdge(0)] = oldtri;
+ trian->tri.push_back(t11);
+
+ ed2update.push_back(t10->getEdge(0));
+ ed2update.push_back(t11->getEdge(0));
+
+ }
+ else{
+ std::list<GEdge*> gGEdges = trian->my_GEdges;
+ std::list<GEdge*>::iterator oe = gGEdges.begin();
+ while(oe != gGEdges.end() && ed.getVertex(0)->onWhat() != *oe && ed.getVertex(0)->onWhat() != (*oe)->getBeginVertex())
+ ++oe;
+ GEdge* ged = *oe;
+ std::vector<MLine*> mlines = ged->lines;
+ std::vector<MLine*> conca;
+ for (unsigned int i=0; i<mlines.size(); i++){
+ if ((mlines[i]->getVertex(0)==ed.getVertex(0)&&mlines[i]->getVertex(1)==ed.getVertex(1))||
+ (mlines[i]->getVertex(1)==ed.getVertex(0)&&mlines[i]->getVertex(0)==ed.getVertex(1))){
+ conca.push_back(new MLine(mlines[i]->getVertex(0),mv));
+ conca.push_back(new MLine(mv,mlines[i]->getVertex(1)));// delete MLine ?
+ }
+ else conca.push_back(mlines[i]);
+ }
+ discreteEdge* de = new discreteEdge(this->model(),NEWLINE(),ged->getBeginVertex(),ged->getEndVertex());
+ setupDiscreteEdge(de,conca,NULL);
+ trian->my_GEdges.remove(ged);
+ ged->lines.clear();
+ ged->mesh_vertices.clear();
+ trian->my_GEdges.push_back(de);
+ trian->borderEdg.erase(ed);
+ trian->borderEdg.insert(t00->getEdge(2));
+ trian->borderEdg.insert(t01->getEdge(1));
+ }
+}
+
+void discreteFace::updateEd2refine(triangulation* trian,std::map<MEdge,double,Less_Edge>&ed2angle,std::vector<MEdge>&ed2update,std::map<double,std::vector<MEdge> >&ed2refine){
+
+ for(unsigned int i=0; i<ed2update.size(); i++){
+ MEdge current = ed2update[i];
+ std::vector<int> intri = trian->ed2tri[current];
+ if(intri.size()>1 && ed2angle.find(current) == ed2angle.end() ){
+ ed2angle[current] = 0.;
+ for(unsigned int j=0; j<intri.size(); j++){
+ MElement* tri = trian->tri[intri[j]];
+ int num = edgeLocalNum(tri,current);
+ MEdge ed1 = tri->getEdge(num+1>2?num-2:num+1);
+ MEdge ed2 = tri->getEdge(num+2>2?num-1:num+2);
+ SVector3 v1(ed1.getVertex(0)->x()-ed1.getVertex(1)->x(),
+ ed1.getVertex(0)->y()-ed1.getVertex(1)->y(),
+ ed1.getVertex(0)->z()-ed1.getVertex(1)->z());
+ SVector3 v2(ed2.getVertex(1)->x()-ed2.getVertex(0)->x(),
+ ed2.getVertex(1)->y()-ed2.getVertex(0)->y(),
+ ed2.getVertex(1)->z()-ed2.getVertex(0)->z());
+ double dotp = dot(v1,v2);
+ double lp = ed1.length()*ed2.length();
+ if(ed2angle[current]==0.)
+ ed2angle[current] = std::acos(dotp/lp);
+ else{
+ ed2angle[current] += std::acos(dotp/lp);
+ if(ed2angle[current] > M_PI)
+ ed2refine[current.length()].push_back(current);
+ else
+ ed2angle.erase(current);
+ }
+ }//end for j
+ }
+ }// end for i
+
+}
+
+
+void discreteFace::updateEd2refineBis(triangulation* trian,std::map<MEdge,double,Less_Edge>&ed2angle,std::vector<MEdge>&ed2update,std::map<double,std::vector<MEdge> >&ed2refine){
+
+ for(unsigned int i=0; i<ed2update.size(); i++){
+ MEdge current = ed2update[i];
+ std::vector<int> intri = trian->ed2tri[current];
+ for(unsigned int j=0; j<intri.size(); j++){
+ if(ed2angle.find(current) == ed2angle.end() ){
+ MElement* tri = trian->tri[intri[j]];
+ int num = edgeLocalNum(tri,current);
+ MEdge ed1 = tri->getEdge(num+1>2?num-2:num+1);
+ MEdge ed2 = tri->getEdge(num+2>2?num-1:num+2);
+ SVector3 v1(ed1.getVertex(0)->x()-ed1.getVertex(1)->x(),
+ ed1.getVertex(0)->y()-ed1.getVertex(1)->y(),
+ ed1.getVertex(0)->z()-ed1.getVertex(1)->z());
+ SVector3 v2(ed2.getVertex(1)->x()-ed2.getVertex(0)->x(),
+ ed2.getVertex(1)->y()-ed2.getVertex(0)->y(),
+ ed2.getVertex(1)->z()-ed2.getVertex(0)->z());
+ double dotp = dot(v1,v2);
+ double lp = ed1.length()*ed2.length();
+ ed2angle[current] = std::acos(dotp/lp);
+ if (ed2angle[current]>M_PI/2.)
+ ed2refine[current.length()].push_back(current);
+ else
+ ed2angle.erase(current);
+ }//end for j
+ }
+ }// end for i
+
+}
+
void discreteFace::addTriangle(triangulation* trian, MTriangle* t)
{
#if defined(HAVE_SOLVER) && defined(HAVE_ANN)
- // #mark quid borders ?
for(int i=0; i<3; i++){
MEdge ed = t->getEdge(i);
- trian->ed2tri[ed].push_back(trian->tri.size());
+ int n = trian->tri.size();
+ trian->fillingHoles.insert(n);
+ trian->ed2tri[ed].push_back(n);
}
trian->tri.push_back(t);
#endif
diff --git a/Geo/discreteFace.h b/Geo/discreteFace.h
index 1f6b0f61af..2515f1fbbc 100644
--- a/Geo/discreteFace.h
+++ b/Geo/discreteFace.h
@@ -37,6 +37,13 @@ class discreteFace : public GFace {
void updateTopology(std::vector<triangulation*>&);
void split(triangulation*,std::vector<triangulation*>&,int);
void fillHoles(triangulation*);
+ void bisectionEdg(triangulation*);
+ bool checkEdges(triangulation*,std::map<double,std::vector<MEdge> >&,std::map<MEdge,double,Less_Edge>&);
+ bool checkEdgesBis(triangulation*,std::map<double,std::vector<MEdge> >&,std::map<MEdge,double,Less_Edge>&);
+ MEdge lepp(triangulation*,MEdge);
+ void refineTriangles(triangulation*,MEdge,std::vector<MEdge>&);
+ void updateEd2refine(triangulation*,std::map<MEdge,double,Less_Edge>&,std::vector<MEdge>&,std::map<double,std::vector<MEdge> >&);
+ void updateEd2refineBis(triangulation*,std::map<MEdge,double,Less_Edge>&,std::vector<MEdge>&,std::map<double,std::vector<MEdge> >&);
void addTriangle(triangulation*,MTriangle*);
GPoint point(double par1, double par2) const;
SPoint2 parFromPoint(const SPoint3 &p, bool onSurface=true) const;
--
GitLab