diff --git a/Geo/discreteDiskFace.cpp b/Geo/discreteDiskFace.cpp
index 3f205ee6b46fe37c9358f4a95e4bcbc279c8f9d5..b7724e6bea6804d6444d876541dacbdc4dfb7d8d 100644
--- a/Geo/discreteDiskFace.cpp
+++ b/Geo/discreteDiskFace.cpp
@@ -3,7 +3,7 @@
// See the LICENSE.txt file for license information. Please report all
// bugs and problems to the public mailing list <gmsh@geuz.org>.
-#include <queue> // #mark
+#include <queue>
#include <stdlib.h>
#include "GmshConfig.h"
@@ -19,14 +19,12 @@
#include "Context.h"
#include "OS.h"
#include "ANN/ANN.h"
-#include "linearSystemCSR.h"
-#include "dofManager.h"
#include "laplaceTerm.h"
#include "convexLaplaceTerm.h"
-#include "convexCombinationTerm.h" // #FIXME
+#include "convexCombinationTerm.h"
#if defined(HAVE_MESH)
-#include "qualityMeasuresJacobian.h" // #temporary?
+#include "qualityMeasuresJacobian.h"
#endif
#if defined(HAVE_OPTHOM)
@@ -286,20 +284,18 @@ discreteDiskFace::discreteDiskFace(GFace *gf, triangulation* diskTriangulation,
_totLength = geoTriangulation->bord.rbegin()->first;
_U0 = geoTriangulation->bord.rbegin()->second;
orderVertices(_totLength, _U0, _coords);
- parametrize(false);
+ parametrize();
buildOct(CAD);
-
-
- if (!checkOrientationUV()){
- Msg::Info("discreteDiskFace:: parametrization is not one-to-one; fixing "
- "the discrete system.");
- parametrize(true);
+ printParamMesh();
+ if (_order > 1 && !checkOrientationUV()){
+ Msg::Info("discreteDiskFace:: parametrization is not one-to-one; the parametric mesh is going to be corrected.");
+ optimize();
buildOct(CAD);
+ printParamMesh();
}
putOnView(true,false);
- printParamMesh();
- if(!checkOrientationUV()) Msg::Fatal("discreteDiskFace:: failing to fix the discrete system");
- // else Msg::Info("Parameterization done :-)");
+ if(!checkOrientationUV()) Msg::Fatal("discreteDiskFace:: failing to compute a one-to-one mapping");
+ else Msg::Info("Parameterization done :-)");
}
/*end BUILDER*/
@@ -320,36 +316,32 @@ void discreteDiskFace::buildOct(std::vector<GFace*> *CAD) const
const int maxElePerBucket = 15;
oct = Octree_Create(maxElePerBucket, origin, ssize, discreteDiskFaceBB,
discreteDiskFaceCentroid, discreteDiskFaceInEle);
-
+
_ddft = new discreteDiskFaceTriangle[discrete_triangles.size()/*-geoTriangulation->fillingHoles.size()*/];
int c = 0;
for(unsigned int i = 0; i < discrete_triangles.size(); ++i){
// 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)];
+ 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);
}
-bool discreteDiskFace::parametrize(bool one2one) const
+bool discreteDiskFace::parametrize()// const
{ // #improveme
- if (one2one)
- Msg::Info("Parametrizing surface %d with 'one-to-one map'", tag());
- else
- Msg::Info("Parametrizing surface %d with 'harmonic map'", tag());
- linearSystem<double> *lsys_u;
- linearSystem<double> *lsys_v;
+ linearSystemCSRTaucs<double> *lsys_u;
+ linearSystemCSRTaucs<double> *lsys_v;
lsys_u = new linearSystemCSRTaucs<double>; // taucs :-)
lsys_v = new linearSystemCSRTaucs<double>; // taucs :-)
@@ -357,20 +349,21 @@ bool discreteDiskFace::parametrize(bool one2one) const
dofManager<double> myAssemblerU(lsys_u); // hashing
dofManager<double> myAssemblerV(lsys_v);
- for(size_t i = 0; i < _U0.size(); i++){
- MVertex *v = _U0[i];
- const double theta = 2 * M_PI * _coords[i];
- 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))])));
+ if(_order>1){
+ for(size_t i = 0; i < _U0.size(); i++){
+ MVertex *v = _U0[i];
+ const double theta = 2 * M_PI * _coords[i];
+ 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){
MElement *t = discrete_triangles[i];
for(int j=0; j<t->getNumVertices(); j++){
@@ -386,28 +379,21 @@ bool discreteDiskFace::parametrize(bool one2one) const
double t1 = Cpu();
simpleFunction<double> ONE(1.0);
+ laplaceTerm mappingU(0, 1, &ONE);
+ laplaceTerm mappingV(0, 1, &ONE);
- if (one2one){
- convexLaplaceTerm mappingU(0, 1, &ONE);
- convexLaplaceTerm mappingV(0, 1, &ONE);
-
- for(unsigned int i = 0; i < discrete_triangles.size(); ++i){
- SElement se(discrete_triangles[i]);
- mappingU.addToMatrix(myAssemblerU, &se);
- mappingV.addToMatrix(myAssemblerV, &se);
- }
+ for(unsigned int i = 0; i < discrete_triangles.size(); ++i){
+ SElement se(discrete_triangles[i]);
+ mappingU.addToMatrix(myAssemblerU, &se);
+ mappingV.addToMatrix(myAssemblerV, &se);
}
- else{
- laplaceTerm mappingU(0, 1, &ONE);
- laplaceTerm mappingV(0, 1, &ONE);
- for(unsigned int i = 0; i < discrete_triangles.size(); ++i){
- SElement se(discrete_triangles[i]);
- mappingU.addToMatrix(myAssemblerU, &se);
- mappingV.addToMatrix(myAssemblerV, &se);
- }
- }
+ if(_order==1){
+ checklsys(lsys_u,&myAssemblerU,1);
+ checklsys(lsys_v,&myAssemblerV,0);
+ }
+
double t2 = Cpu();
Msg::Debug("Assembly done in %8.3f seconds", t2 - t1);
lsys_u->systemSolve();
@@ -439,6 +425,91 @@ bool discreteDiskFace::parametrize(bool one2one) const
return true;
}
+void discreteDiskFace::checklsys(linearSystemCSRTaucs<double>* lsys,dofManager<double>* myAssembler,int U)
+{
+
+ INDEX_TYPE *ai;
+ INDEX_TYPE *jptr;
+ double *a;
+ //double b;
+ lsys->getMatrix(jptr,ai,a);
+ int NbUnk = lsys->getNbUnk();
+ // int NNZ = lsys->getNNZ();
+ double tr = 0.;
+ //printf("------------------------%d\n######################%d\n______________________\n",NNZ,NbUnk);
+ for(int i=0; i<NbUnk; i++){
+ tr += a[jptr[i]];
+ /*lsys->getFromRightHandSide(i,b);
+ printf("%d: %d->%d\n",i,jptr[i],jptr[i+1]);
+ for(int j=jptr[i]; j<jptr[i+1]; j++)
+ printf("%d: %f \t",ai[j],a[j]);
+ printf("=%f \n",b);*/
+ }
+ tr /= NbUnk;
+
+ double alpha = .1;
+ for(int i=0; i<NbUnk; i++)
+ for(int j=jptr[i]+1; j<jptr[i+1]; j++)
+ if(a[j] > -alpha*tr){
+ double delta = a[j] + alpha*tr;
+ a[jptr[i]] += delta;
+ a[jptr[ai[j]]] += delta;
+ a[j] -= delta;
+ }
+ /*
+ printf("\n");
+ for(int i=0; i<NbUnk; i++){
+ lsys->getFromRightHandSide(i,b);
+ printf("%d: %d->%d\n",i,jptr[i],jptr[i+1]);
+ for(int j=jptr[i]; j<jptr[i+1]; j++)
+ printf("%d: %f \t",ai[j],a[j]);
+ printf("=%f \n",b);
+ }
+ */
+ std::map<int,int> dirichlet;
+ for(unsigned int i=0; i<_U0.size(); i++)
+ dirichlet[myAssembler->getDofNumber(Dof(_U0[i]->getNum(),10000))] = i;
+
+
+ for(int i=0; i<NbUnk; i++){
+ if(dirichlet.find(i)!=dirichlet.end()){
+ lsys->addToMatrix(i,i,-a[jptr[i]]+1);
+ double theta = 2 * M_PI * _coords[dirichlet[i]];
+ double val = U*cos(theta) + (1-U)*sin(theta);
+ lsys->addToRightHandSide(i,val);
+ for(int j=jptr[i]+1; j<jptr[i+1]; j++){
+ if(dirichlet.find(ai[j])==dirichlet.end())
+ lsys->addToRightHandSide(ai[j],-a[j]*val);
+ a[j] = 0;
+ }
+ }
+ else{
+ for(int j=jptr[i]+1; j<jptr[i+1]; j++){
+ if(dirichlet.find(ai[j])!=dirichlet.end()){
+ double theta = 2 * M_PI * _coords[dirichlet[ai[j]]];
+ double val = U*cos(theta) + (1-U)*sin(theta);
+ lsys->addToRightHandSide(i,-a[j]*val);
+ a[j] = 0;
+ }
+ }
+ }//end else
+ }//end for i
+ /*
+ lsys->getMatrix(jptr,ai,a);
+ NbUnk = lsys->getNbUnk();
+ NNZ = lsys->getNNZ();
+ printf("--#--#--#--\n");
+ for(int i=0; i<NbUnk; i++){
+ lsys->getFromRightHandSide(i,b);
+ printf("%d: %d->%d\n",i,jptr[i],jptr[i+1]);
+ for(int j=jptr[i]; j<jptr[i+1]; j++)
+ printf("%d: %f \t",ai[j],a[j]);
+ printf("=%f \n",b);
+ }
+ */
+
+}
+
void discreteDiskFace::getTriangleUV(const double u,const double v,
discreteDiskFaceTriangle **mt,
double &_xi, double &_eta)const{
@@ -479,7 +550,7 @@ bool discreteDiskFace::checkOrientationUV()
else nbP++;
}
if (nbP*nbM){
- Msg::Error("Map %d of the atlas : Triangles have different orientations (%d + / %d -)",tag(),nbP,nbM);
+ Msg::Info("Map %d of the atlas : Triangles have different orientations (%d + / %d -)",tag(),nbP,nbM);
return false;
}
return true;
@@ -508,17 +579,20 @@ bool discreteDiskFace::checkOrientationUV()
}
void discreteDiskFace::optimize()
-{
+{ // #improve . . .
#if defined(HAVE_OPTHOM)
// parameters for mesh optimization
// -- high order
OptHomParameters optParams;
optParams.dim = 2;
- optParams.optPassMax = 10;
+ optParams.optPassMax = 100;
optParams.TMAX = 300;
optParams.fixBndNodes = true;
- optParams.strategy = 1;
+ optParams.optPrimSurfMesh = true;
+ optParams.BARRIER_MIN = 1e-3;
+ optParams.BARRIER_MAX = 1e3;
+ optParams.strategy = 0;
// -- linear
MeshQualOptParameters opt;
opt.dim = 2;
@@ -544,32 +618,47 @@ void discreteDiskFace::optimize()
paramTriangles[i] = new MTriangle6(mv);
}
// -- generation of the parametric topology #mark what about the GFace for the GModel ?
- std::map<int,std::vector<MElement*> > e2e;
+ std::map<int,std::vector<MElement*> > e2e;// entity to element
e2e[0] = paramTriangles;
std::vector<int> v;
v.push_back(0);
- std::map<int,std::vector<int> > e2p;
+ std::map<int,std::vector<int> > e2p;// entity to physical
e2p[0] = v;
+ std::set<MVertex*> todelete;
GModel* paramDisk = GModel::createGModel(e2e,e2p);
- discreteVertex dv(paramDisk,NEWPOINT());
- dv.mesh_vertices.push_back(sp2mv[coordinates[_U0[0]]]);
- sp2mv[coordinates[_U0[0]]]->setEntity(&dv);
- paramDisk->add(&dv);
- discreteEdge de(paramDisk,NEWLINE(),&dv,&dv);
- paramDisk->add(&de);
+ // ---- discrete vertex
+ std::set<GFace*, GEntityLessThan>::iterator it = paramDisk->firstFace();
+ GFace *dgf = *it;
+ discreteVertex *dv = new discreteVertex(paramDisk,paramDisk->getMaxElementaryNumber(0)+1);
+ sp2mv[coordinates[_U0[0]]]->setEntity(dv);
+ dv->mesh_vertices.push_back(sp2mv[coordinates[_U0[0]]]);
+ todelete.insert(sp2mv[coordinates[_U0[0]]]);
+ paramDisk->add(dv);
+ // ---- discrete edge
+ discreteEdge *de = new discreteEdge(paramDisk,paramDisk->getMaxElementaryNumber(1)+1,dv,dv);
+ paramDisk->add(de);
+ std::vector<MLine*> lines;
for(unsigned int i=1; i<_U0.size(); i++){
- sp2mv[coordinates[_U0[i]]]->setEntity(&de);
- de.mesh_vertices.push_back(sp2mv[coordinates[_U0[i]]]);
- de.lines.push_back(new MLine(sp2mv[coordinates[_U0[i-1]]],sp2mv[coordinates[_U0[i]]]));
+ sp2mv[coordinates[_U0[i]]]->setEntity(de);
+ de->mesh_vertices.push_back(sp2mv[coordinates[_U0[i]]]);
+ todelete.insert(sp2mv[coordinates[_U0[i]]]);
+ lines.push_back(new MLine(sp2mv[coordinates[_U0[i-1]]],sp2mv[coordinates[_U0[i]]]));
}
- de.createGeometry();
-
+ lines.push_back(new MLine(sp2mv[coordinates[_U0[_U0.size()-1]]],sp2mv[coordinates[_U0[0]]]));
+ de->setTopo(lines);
+ de->createGeometry();// !!!! setTopo ... MLine's
+
+
+
// optimization
if(_order >1)
HighOrderMeshOptimizer(paramDisk, optParams);
else
MeshQualityOptimizer(paramDisk,opt);
+
+
+
// update the parametrization
paramTriangles = e2e[0];
for(unsigned int i=0; i< paramTriangles.size(); i++){
@@ -583,8 +672,19 @@ void discreteDiskFace::optimize()
}
}
+ // update GFace's mesh_vertices
+ std::vector<MVertex*> newMV;
+ for(unsigned int imv=0; imv<dgf->mesh_vertices.size(); imv++){
+ MVertex* current = dgf->mesh_vertices[imv];
+ std::set<MVertex*>::iterator itmv = todelete.find(current);
+ if (itmv==todelete.end()) newMV.push_back(current);
+ }
+ dgf->mesh_vertices.clear();
+ dgf->mesh_vertices = newMV;
+
// cleaning
delete paramDisk;
+
#endif
}
@@ -738,7 +838,7 @@ void discreteDiskFace::secondDer(const SPoint2 ¶m,
}
void discreteDiskFace::putOnView(bool Xu, bool Ux)
-{
+{// #improveme using built-in methods
char mybuffer [64];
@@ -979,35 +1079,6 @@ GPoint discreteDiskFace::intersectionWithCircle(const SVector3 &n1, const SVecto
return pp;
}
-void discreteDiskFace::printParamMesh()
-{
- std::map<MVertex*,int> mv2int;
- char buffer[256];
- sprintf(buffer,"param_mesh%d.msh",tag());
- FILE* pmesh = fopen(buffer,"w");
- int count = 1;
- fprintf(pmesh,"$MeshFormat\n2.2 0 8\n$EndMeshFormat\n$Nodes\n%u\n",(unsigned int)allNodes.size());
- for(std::set<MVertex*>::iterator it = allNodes.begin(); it!=allNodes.end(); ++it){
- fprintf(pmesh,"%d %f %f 0\n",count,(coordinates[(*it)]).x(),(coordinates[(*it)]).y());
- mv2int[*it] = count;
- count++;
- }
- fprintf(pmesh,"$EndNodes\n$Elements\n%u\n",(unsigned int)discrete_triangles.size());
- for(unsigned int i=0; i<discrete_triangles.size(); i++){
- MElement* tri = discrete_triangles[i];
- int p;
- _order == 1 ? p=2 : p=9;
- fprintf(pmesh,"%d %d 2 0 0",i+1,p);
- for(int j=0; j<_N; j++){
- MVertex* mv = tri->getVertex(j);
- fprintf(pmesh," %d",mv2int[mv]);
- }
- fprintf(pmesh,"\n");
- }
- fprintf(pmesh,"$EndElements\n");
- fclose(pmesh);
-}
-
// computes some kind of maximal distance in a mesh
static void addTo (std::map<MVertex*, std::vector<MElement*> > &v2t, MVertex *v, MElement *t)
@@ -1036,7 +1107,7 @@ static MEdge getEdge (MElement *t, MVertex *v)
return MEdge();
}
-/* warning
+/* #warning
static double computeDistanceLinePoint (MVertex *v1, MVertex *v2, MVertex *v){
SVector3 U = v2->point() - v1->point();
@@ -1196,4 +1267,71 @@ double triangulation::geodesicDistance ()
return CLOSEST;
}
+void discreteDiskFace::printAtlasMesh()
+{
+#if defined(HAVE_SOLVER) && defined(HAVE_ANN)
+ std::map<MVertex*,int> mv2int;
+ char buffer[256];
+ sprintf(buffer,"atlas_mesh%d.msh",initialTriangulation->idNum);
+ FILE* pmesh = Fopen(buffer,"w");
+
+ std::set<MVertex*> meshvertices;
+
+ for(unsigned int i=0; i<initialTriangulation->tri.size(); ++i){
+ MElement* tri = initialTriangulation->tri[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",(unsigned int)initialTriangulation->tri.size());
+ for(unsigned int i=0; i<initialTriangulation->tri.size(); i++){
+ MElement* tri = initialTriangulation->tri[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);
+#endif
+}
+
+void discreteDiskFace::printParamMesh()
+{
+ std::map<MVertex*,int> mv2int;
+ char buffer[256];
+ sprintf(buffer,"param_mesh%d.msh",tag());
+ FILE* pmesh = fopen(buffer,"w");
+ int count = 1;
+ fprintf(pmesh,"$MeshFormat\n2.2 0 8\n$EndMeshFormat\n$Nodes\n%u\n",(unsigned int)allNodes.size());
+ for(std::set<MVertex*>::iterator it = allNodes.begin(); it!=allNodes.end(); ++it){
+ fprintf(pmesh,"%d %f %f 0\n",count,(coordinates[(*it)]).x(),(coordinates[(*it)]).y());
+ mv2int[*it] = count;
+ count++;
+ }
+ fprintf(pmesh,"$EndNodes\n$Elements\n%u\n",(unsigned int)discrete_triangles.size());
+ for(unsigned int i=0; i<discrete_triangles.size(); i++){
+ MElement* tri = discrete_triangles[i];
+ int p;
+ _order == 1 ? p=2 : p=9;
+ fprintf(pmesh,"%d %d 2 0 0",i+1,p);
+ for(int j=0; j<_N; j++){
+ MVertex* mv = tri->getVertex(j);
+ fprintf(pmesh," %d",mv2int[mv]);
+ }
+ fprintf(pmesh,"\n");
+ }
+ fprintf(pmesh,"$EndElements\n");
+ fclose(pmesh);
+}
+
#endif
diff --git a/Geo/discreteDiskFace.h b/Geo/discreteDiskFace.h
index 170d87df7f9a57107a9b284351d330edd78339f0..b5aed56fa930f0c7791b5ff736f73b2056aeda80 100644
--- a/Geo/discreteDiskFace.h
+++ b/Geo/discreteDiskFace.h
@@ -23,6 +23,8 @@
#include "PView.h"
#include "robustPredicates.h"
#include "MLine.h"
+#include "linearSystemCSR.h"
+#include "dofManager.h"
inline int nodeLocalNum(MElement* e, MVertex* v)
{
@@ -145,7 +147,7 @@ class triangulation {
std::map<MVertex*,std::vector<MVertex*> >::iterator im = firstNode2Edge.find(first);
if (im != firstNode2Edge.end()){
- Msg::Info("Supposed seam point for discreteFace %d", gf->tag());
+ Msg::Info("Assumed seam point for discreteFace %d", gf->tag());
seamPoint = true;
break;
}
@@ -198,7 +200,7 @@ class triangulation {
fprintf(f,"};\n");
fclose(f);
}
-
+
triangulation() : gf(0) {}
triangulation(int id, std::vector<MElement*> input, GFace* gface)
: idNum(id), tri(input), gf(gface) { assign(); }
@@ -217,11 +219,11 @@ class discreteDiskFaceTriangle {
class discreteDiskFace : public GFace {
GFace *_parent;
void buildOct(std::vector<GFace*> *CAD = NULL) const;
- bool parametrize(bool one2one = false) const;
+ bool parametrize();// const;
+ void checklsys(linearSystemCSRTaucs<double>*,dofManager<double>*,int);
void putOnView(bool,bool);
bool checkOrientationUV();
void optimize();
- void printParamMesh();
public:
discreteDiskFace(GFace *parent, triangulation* diskTriangulation,
@@ -241,14 +243,19 @@ class discreteDiskFace : public GFace {
GPoint intersectionWithCircle(const SVector3 &n1, const SVector3 &n2,
const SVector3 &p, const double &d,
double uv[2]) const;
+ void printAtlasMesh () ;
+ void printParamMesh () ;
+
std::vector<MElement*> discrete_triangles;
+
+
+
protected:
// a copy of the mesh that should not be destroyed
triangulation* initialTriangulation;
triangulation* geoTriangulation;// parametrized triangulation
-
std::vector<MVertex*> discrete_vertices;
int _order;
diff --git a/Geo/discreteFace.cpp b/Geo/discreteFace.cpp
index 0f94d4cd450dcd89c93b2a577e9fcf4912ae9513..ea1c672c1e31e972dd201ba14370f8b02f99728d 100644
--- a/Geo/discreteFace.cpp
+++ b/Geo/discreteFace.cpp
@@ -14,6 +14,8 @@
#include <stack>
#include <queue>
+#include "MPoint.h"
+
#if defined(HAVE_METIS)
extern "C" {
#include <metis.h>
@@ -109,8 +111,9 @@ void discreteFace::secondDer(const SPoint2 ¶m,
void discreteFace::createGeometry()
{
- checkAndFixOrientation();
+ checkAndFixOrientation();
+
#if defined(HAVE_SOLVER) && defined(HAVE_ANN)
int order = 1;
int nPart = 2;
@@ -121,19 +124,16 @@ void discreteFace::createGeometry()
int id=1;
std::stack<triangulation*> toSplit;
std::vector<triangulation*> toParam;
+
std::vector<MElement*> tem(triangles.begin(),triangles.end());
-
triangulation* init = new triangulation(-1, tem,this);
toSplit.push(init);
if((toSplit.top())->genus()!=0 || (toSplit.top())->aspectRatio() > eta ||
- (toSplit.top())->seamPoint){
-
+ (toSplit.top())->seamPoint){
while( !toSplit.empty()){
std::vector<triangulation*> part;
triangulation* tosplit = toSplit.top();
toSplit.pop();
- // printf("genus %d ar %12.5E\n",tosplit->genus(), tosplit->aspectRatio());
- // getchar();
split(tosplit,part,nPart);
delete tosplit;
@@ -143,7 +143,6 @@ void discreteFace::createGeometry()
else{
toParam.push_back(part[i]);
part[i]->idNum=id++;
- // printf("part %d is OK\n", part[i]->idNum);
}
}// end for i
}// !.empty()
@@ -152,24 +151,23 @@ void discreteFace::createGeometry()
toParam.push_back(toSplit.top());
toSplit.top()->idNum=id++;
}
- double t1 = Cpu();
- printf("split done\n");
updateTopology(toParam);
- double t2 = Cpu();
- printf("topology done in %8.3f sec\n",t2-t1);
-
+
for(unsigned int i=0; i<toParam.size(); i++){
- printf("MAP(%d) : aspect ratio = %12.5E\n",toParam[i]->idNum,toParam[i]->aspectRatio());
+ /*
char name[256];
sprintf(name,"map%d.pos",i);
toParam[i]->print(name,i);
+ */
fillHoles(toParam[i]);
- sprintf(name,"mapFilled%d.pos",i);
- toParam[i]->print(name, toParam[i]->idNum);
+ //sprintf(name,"mapFilled%d.pos",i);
+ //toParam[i]->print(name, toParam[i]->idNum);
}
- 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->printAtlasMesh();
df->replaceEdges(toParam[i]->my_GEdges);
_atlas.push_back(df);
}
@@ -224,96 +222,26 @@ 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++){
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);
+ _atlas[i]->mesh(verbose);/*
+ const char *name = "atlas%d";
+ char filename[256];
+ sprintf(filename,name,i);
+ _atlas[i]->printPhysicalMesh(filename);*/
}
-
+
gatherMeshes();
meshStatistics.status = GFace::DONE;
#endif
}
-void discreteFace::printAtlasMesh(std::vector<MElement*> elm, int I)
-{
-
- std::map<MVertex*,int> mv2int;
- char buffer[256];
- 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",(unsigned int)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)
-{
-#if defined(HAVE_SOLVER) && defined(HAVE_ANN)
- 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",(unsigned int)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);
-#endif
-}
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)
@@ -393,10 +321,12 @@ void discreteFace::checkAndFixOrientation()
} // end while
}
+
void discreteFace::setupDiscreteVertex(GVertex*dv,MVertex*mv,std::set<MVertex*>*trash){
- dv->mesh_vertices.push_back(mv);
mv->setEntity(dv);
+ dv->mesh_vertices.push_back(mv);
this->model()->add(dv);
+ dv->points.push_back(new MPoint(dv->mesh_vertices.back()));
if (trash) trash->insert(mv);
}
@@ -404,15 +334,16 @@ void discreteFace::setupDiscreteEdge(discreteEdge*de,std::vector<MLine*>mlines,
std::set<MVertex*>*trash)
{
this->model()->add(de);// new GEdge
- de->lines = mlines;// associated MLine's
+ de->setTopo(mlines);// associated MLine's
for(unsigned int i=1; i<mlines.size(); i++){//not the first vertex of the GEdge (neither the last one)
mlines[i]->getVertex(0)->setEntity(de);
de->mesh_vertices.push_back(mlines[i]->getVertex(0));
if(trash) trash->insert(mlines[i]->getVertex(0));
}
de->createGeometry();
- if(de->getBeginVertex()->mesh_vertices[0]!=de->lines[0]->getVertex(0))// small hack ... #mark
- de->reverse();
+ de->getBeginVertex()->addEdge(de);
+ de->getEndVertex()->addEdge(de);
+
}
// split old GEdge's
@@ -420,9 +351,9 @@ void discreteFace::splitDiscreteEdge(GEdge *de , GVertex *gv, discreteEdge* newE
{
MVertex *vend = de->getEndVertex()->mesh_vertices[0];
- newE[0] = new discreteEdge (de->model(),NEWLINE(),de->getBeginVertex(),gv);
- newE[1] = new discreteEdge (de->model(),NEWLINE(),gv, de->getEndVertex());
-
+ newE[0] = new discreteEdge (de->model(),model()->getMaxElementaryNumber(1) + 1,de->getBeginVertex(),gv);
+ newE[1] = new discreteEdge (de->model(),model()->getMaxElementaryNumber(1) + 1,gv, de->getEndVertex());
+
int current = 0;
std::vector<MLine*> mlines;
// assumption: the MLine's are sorted
@@ -435,6 +366,8 @@ void discreteFace::splitDiscreteEdge(GEdge *de , GVertex *gv, discreteEdge* newE
current++;
}
}
+ de->getBeginVertex()->delEdge(de);
+ de->getEndVertex()->delEdge(de);
de->mesh_vertices.clear();
de->lines.clear();
@@ -569,7 +502,7 @@ void discreteFace::updateTopology(std::vector<triangulation*>&partition)
std::set<MVertex*> todelete; // vertices that do not belong to the GFace anymore
std::set<GEdge*> gGEdges(l_edges.begin(),l_edges.end());// initial GEdges of the GFace (to be updated)
- for(int i=0; i<nPartitions; i++){// each part is going ot be compared with the other ones
+ for(int i=0; i<nPartitions; i++){// each part is going to be compared with the other ones
const 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 partitions, with ii > i since ii < i have already been compared
std::map<MVertex*,MLine*> v02edg;//first vertex of the MLine
@@ -606,25 +539,25 @@ void discreteFace::updateTopology(std::vector<triangulation*>&partition)
mvt.resize(2);
mvt[0] = *itf;
mvt[1] = cv0;
-
// creation of the GVertex, for new nonloop GEdge's
for(int imvt=0; imvt<2; imvt++){
std::set<GEdge*>::iterator oe=gGEdges.begin();// find the old GEdge that has the current new GVertex
while(oe !=gGEdges.end() && mvt[imvt]->onWhat() != *oe && mvt[imvt]->onWhat() != (*oe)->getBeginVertex() && mvt[imvt]->onWhat() != (*oe)->getEndVertex())
- ++oe;
+ ++oe;
+
if (oe == gGEdges.end()){// not on an existing GEdge: new internal GVertex
- v[imvt] = new discreteVertex (this->model(),NEWPOINT());
+ v[imvt] = new discreteVertex (this->model(),model()->getMaxElementaryNumber(0) + 1);
setupDiscreteVertex(v[imvt],mvt[imvt],&todelete);
}
else{// on an existing GEdge
GEdge* onit = *oe;// the new GVertex can already exist; if it is the case, there is no need to create a new one
if(mvt[imvt] == onit->getBeginVertex()->mesh_vertices[0])
- v[imvt] = onit->getBeginVertex();
+ v[imvt] = onit->getBeginVertex();
else if (mvt[imvt] == onit->getEndVertex()->mesh_vertices[0])
- v[imvt] = onit->getEndVertex();
+ v[imvt] = onit->getEndVertex();
else{
- v[imvt] = new discreteVertex (this->model(),NEWPOINT());
+ v[imvt] = new discreteVertex (this->model(),model()->getMaxElementaryNumber(0) + 1);
setupDiscreteVertex(v[imvt],mvt[imvt],NULL);
discreteEdge* de[2];
splitDiscreteEdge(onit,v[imvt],de);
@@ -635,11 +568,9 @@ void discreteFace::updateTopology(std::vector<triangulation*>&partition)
}// end else oe==end()
}// end imvt
// the new GEdge can be created with its GVertex
- discreteEdge* internalE = new discreteEdge (this->model(),NEWLINE(),v[0],v[1]);
- // printf("creating new discrete edge %d with %d lines (%p %p)\n", internalE->tag(), myLines.size(),v[0],v[1]);
+ discreteEdge* internalE = new discreteEdge (this->model(),model()->getMaxElementaryNumber(1) + 1,v[0],v[1]);
setupDiscreteEdge(internalE,myLines,&todelete);
- partition[i]->my_GEdges.push_back(internalE);
- partition[ii]->my_GEdges.push_back(internalE);
+ gGEdges.insert(internalE);
first.erase(itf);// next first vertex of a nonloop GEdge
}//end while first.empty()
@@ -654,17 +585,16 @@ void discreteFace::updateTopology(std::vector<triangulation*>&partition)
v02edg.erase(cv0);// update the container
cv0 = my_MLines.back()->getVertex(1);// next MLine of the loop
}
- v = new discreteVertex (this->model(),NEWPOINT());
+ v = new discreteVertex (this->model(),model()->getMaxElementaryNumber(0) + 1);
setupDiscreteVertex(v,cv0,&todelete);
- discreteEdge* gloop = new discreteEdge (this->model(),NEWLINE(),v,v);
+ discreteEdge* gloop = new discreteEdge (this->model(),model()->getMaxElementaryNumber(1) + 1,v,v);
setupDiscreteEdge(gloop,my_MLines,&todelete);
- partition[i]->my_GEdges.push_back(gloop);
- partition[ii]->my_GEdges.push_back(gloop);
+ gGEdges.insert(gloop);
}// end while v02edg.empty()
}//end for ii
}// end for i
-
-
+
+
// adding old-updated bounding GEdge's to the corresponding partitions
for(std::set<GEdge*>::iterator le=gGEdges.begin(); le!=gGEdges.end(); ++le){
GEdge* ile = *le;
@@ -673,7 +603,7 @@ void discreteFace::updateTopology(std::vector<triangulation*>&partition)
std::set<MEdge,Less_Edge> bordi = partition[i]->borderEdg;
if(bordi.find(edTest)!=bordi.end()){
partition[i]->my_GEdges.push_back(ile);
- break;
+ //break;
}
}
}
@@ -691,7 +621,7 @@ void discreteFace::updateTopology(std::vector<triangulation*>&partition)
}
void discreteFace::fillHoles(triangulation* trian)
-{
+{// #improveme moving the center
#if defined(HAVE_SOLVER) && defined(HAVE_ANN)
std::map<double,std::vector<MVertex*> > bords = trian->bord;
std::map<double,std::vector<MVertex*> >::reverse_iterator it = bords.rbegin();
@@ -749,311 +679,6 @@ void discreteFace::fillHoles(triangulation* trian)
#endif
}
-void discreteFace::bisectionEdg(triangulation* trian)
-{
-#if defined(HAVE_SOLVER) && defined(HAVE_ANN)
- 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
-#endif
-}
-
-bool discreteFace::checkEdges(triangulation* trian,
- std::map<double,std::vector<MEdge> >& ed2refine,
- std::map<MEdge,double,Less_Edge>& ed2angle)
-{
-#if defined(HAVE_SOLVER) && defined(HAVE_ANN)
- 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());
-#else
- return false;
-#endif
-}
-
-
-bool discreteFace::checkEdgesBis(triangulation* trian,
- std::map<double,std::vector<MEdge> >& ed2refine,
- std::map<MEdge,double,Less_Edge>& ed2angle)
-{
-#if defined(HAVE_SOLVER) && defined(HAVE_ANN)
- 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());
-#else
- return false;
-#endif
-}
-
-MEdge discreteFace::lepp(triangulation* trian,MEdge ed)
-{
-#if defined(HAVE_SOLVER) && defined(HAVE_ANN)
- 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;
-#else
- return MEdge();
-#endif
-}
-
-void discreteFace::refineTriangles(triangulation* trian,MEdge ed,
- std::vector<MEdge>&ed2update)
-{
-#if defined(HAVE_SOLVER) && defined(HAVE_ANN)
- 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));
- }
-#endif
-}
-
-void discreteFace::updateEd2refine(triangulation* trian,
- std::map<MEdge,double,Less_Edge>&ed2angle,
- std::vector<MEdge>&ed2update,
- std::map<double,std::vector<MEdge> >&ed2refine)
-{
-#if defined(HAVE_SOLVER) && defined(HAVE_ANN)
- 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
-#endif
-}
-
-void discreteFace::updateEd2refineBis(triangulation* trian,
- std::map<MEdge,double,Less_Edge>&ed2angle,
- std::vector<MEdge>&ed2update,
- std::map<double,std::vector<MEdge> >&ed2refine)
-{
-#if defined(HAVE_SOLVER) && defined(HAVE_ANN)
- 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
-#endif
-}
-
void discreteFace::addTriangle(triangulation* trian, MTriangle* t)
{
#if defined(HAVE_SOLVER) && defined(HAVE_ANN)
diff --git a/Geo/discreteFace.h b/Geo/discreteFace.h
index 2515f1fbbc3ce9326852958597e8223d2a0f4d69..ee911f6e76eefedef6e7da68e3d5a7716c31996a 100644
--- a/Geo/discreteFace.h
+++ b/Geo/discreteFace.h
@@ -31,19 +31,13 @@ class discreteFace : public GFace {
discreteFace(GModel *model, int num);
virtual ~discreteFace() {}
void checkAndFixOrientation();
+ void checkConnectivity(std::vector<std::vector<MElement*> >&);
void setupDiscreteVertex(GVertex*,MVertex*,std::set<MVertex*>*);
void setupDiscreteEdge(discreteEdge*,std::vector<MLine*>,std::set<MVertex*>*);
void splitDiscreteEdge(GEdge*,GVertex*,discreteEdge*[2]);
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;