diff --git a/Geo/GFaceCompound.cpp b/Geo/GFaceCompound.cpp
index 8b9b4ab4e91440ba0f1257ede40f92fc709fafe5..203d8fca12d6d4fd5a365571f791cf9461a02015 100644
--- a/Geo/GFaceCompound.cpp
+++ b/Geo/GFaceCompound.cpp
@@ -44,6 +44,40 @@ static void fixEdgeToValue(GEdge *ed, double value, dofManager<double> &myAssemb
}
}
+//--------------------------------------------------------------
+static int intersection_segments (SPoint3 &p1, SPoint3 &p2,
+ SPoint3 &q1, SPoint3 &q2,
+ double x[2]){
+
+
+ double xp_max = std::max(p1.x(),p2.x());
+ double yp_max = std::max(p1.y(),p2.y());
+ double xq_max = std::max(q1.x(),q2.x());
+ double yq_max = std::max(q1.y(),q2.y());
+
+ double xp_min = std::min(p1.x(),p2.x());
+ double yp_min = std::min(p1.y(),p2.y());
+ double xq_min = std::min(q1.x(),q2.x());
+ double yq_min = std::min(q1.y(),q2.y());
+ if ( yq_min > yp_max || xq_min > xp_max ||
+ yq_max < yp_min || xq_max < xp_min){
+ return 0;
+ }
+ else{
+ double A[2][2];
+ A[0][0] = p2.x()-p1.x();
+ A[0][1] = q1.x()-q2.x();
+ A[1][0] = p2.y()-p1.y();
+ A[1][1] = q1.y()-q2.y();
+ double b[2] = {q1.x()-p1.x(),q1.y()-p1.y()};
+ sys2x2(A,b,x);
+
+ return (x[0] >= 0.0 && x[0] <= 1. &&
+ x[1] >= 0.0 && x[1] <= 1.);
+ }
+
+}
+//--------------------------------------------------------------
static void computeCGKernelPolygon(std::map<MVertex*,SPoint3> &coordinates,
std::vector<MVertex*> &cavV, double &ucg, double &vcg)
{
@@ -242,6 +276,8 @@ void GFaceCompound::fillNeumannBCS() const
double x=0.;
double y=0.;
double z=0.;
+ //EMI- TODO FIND KERNEL OF POLYGON AND PLACE AT CG KERNEL !
+ //IF NO KERNEL -> DO NOT FILL TRIS
for (std::list<GEdge*>::iterator ite = loop.begin(); ite != loop.end(); ite++){
for (unsigned int k= 0; k< (*ite)->getNumMeshElements(); k++){
MVertex *v0 = (*ite)->getMeshElement(k)->getVertex(0);
@@ -261,46 +297,45 @@ void GFaceCompound::fillNeumannBCS() const
MVertex *v0 = (*ite)->getMeshElement(i)->getVertex(0);
MVertex *v1 = (*ite)->getMeshElement(i)->getVertex(1);
- //fillTris.push_back(new MTriangle(v0,v1, c));
+// fillTris.push_back(new MTriangle(v0,v1, c));
-// MVertex *v2 = new MVertex(.5*(v0->x()+c->x()), .5*(v0->y()+c->y()), .5*(v0->z()+c->z()));
-// MVertex *v3 = new MVertex(.5*(v1->x()+c->x()), .5*(v1->y()+c->y()), .5*(v1->z()+c->z()));
-// fillTris.push_back(new MTriangle(v0,v2,v3));
-// fillTris.push_back(new MTriangle(v2,c, v3));
-// fillTris.push_back(new MTriangle(v0,v3, v1);
-
- MVertex *v2 = new MVertex(.66*v0->x()+.33*c->x(), .66*v0->y()+.33*c->y(), .66*v0->z()+.33*c->z());
- MVertex *v3 = new MVertex(.66*v1->x()+.33*c->x(), .66*v1->y()+.33*c->y(), .66*v1->z()+.33*c->z());
- MVertex *v4 = new MVertex(.33*v0->x()+.66*c->x(), .33*v0->y()+.66*c->y(), .33*v0->z()+.66*c->z());
- MVertex *v5 = new MVertex(.33*v1->x()+.66*c->x(), .33*v1->y()+.66*c->y(), .33*v1->z()+.66*c->z());
+ MVertex *v2 = new MVertex(.5*(v0->x()+c->x()), .5*(v0->y()+c->y()), .5*(v0->z()+c->z()));
+ MVertex *v3 = new MVertex(.5*(v1->x()+c->x()), .5*(v1->y()+c->y()), .5*(v1->z()+c->z()));
fillTris.push_back(new MTriangle(v0,v2,v3));
- fillTris.push_back(new MTriangle(v2,v5,v3));
- fillTris.push_back(new MTriangle(v2,v4,v5));
- fillTris.push_back(new MTriangle(v4,c,v5));
- fillTris.push_back(new MTriangle(v0,v3,v1));
+ fillTris.push_back(new MTriangle(v2,c, v3));
+ fillTris.push_back(new MTriangle(v0,v3, v1)) ;
+
+// MVertex *v2 = new MVertex(.66*v0->x()+.33*c->x(), .66*v0->y()+.33*c->y(), .66*v0->z()+.33*c->z());
+// MVertex *v3 = new MVertex(.66*v1->x()+.33*c->x(), .66*v1->y()+.33*c->y(), .66*v1->z()+.33*c->z());
+// MVertex *v4 = new MVertex(.33*v0->x()+.66*c->x(), .33*v0->y()+.66*c->y(), .33*v0->z()+.66*c->z());
+// MVertex *v5 = new MVertex(.33*v1->x()+.66*c->x(), .33*v1->y()+.66*c->y(), .33*v1->z()+.66*c->z());
+// fillTris.push_back(new MTriangle(v0,v2,v3));
+// fillTris.push_back(new MTriangle(v2,v5,v3));
+// fillTris.push_back(new MTriangle(v2,v4,v5));
+// fillTris.push_back(new MTriangle(v4,c,v5));
+// fillTris.push_back(new MTriangle(v0,v3,v1));
}
}
}
}
- //if (fillTris.size() > 0)
- //printf("**** Filling Neuman BCs with %d triangles \n", fillTris.size());
-
-// FILE * ftri = fopen("fillTris.pos","w");
-// fprintf(ftri,"View \"\"{\n");
-// for (std::list<MTriangle*>::iterator it2 = fillTris.begin(); it2 !=fillTris.end(); it2++ ){
-// MTriangle *t = (*it2);
-// fprintf(ftri,"ST(%g,%g,%g,%g,%g,%g,%g,%g,%g){%g,%g,%g};\n",
-// t->getVertex(0)->x(), t->getVertex(0)->y(), t->getVertex(0)->z(),
-// t->getVertex(1)->x(), t->getVertex(1)->y(), t->getVertex(1)->z(),
-// t->getVertex(2)->x(), t->getVertex(2)->y(), t->getVertex(2)->z(),
-// 1., 1., 1.);
-// }
-// fprintf(ftri,"};\n");
-// fclose(ftri);
-
-// exit(1);
+
+ if (fillTris.size() > 0){
+ FILE * ftri = fopen("fillTris.pos","a");
+ fprintf(ftri,"View \"\"{\n");
+ for (std::list<MTriangle*>::iterator it2 = fillTris.begin(); it2 !=fillTris.end(); it2++ ){
+ MTriangle *t = (*it2);
+ fprintf(ftri,"ST(%g,%g,%g,%g,%g,%g,%g,%g,%g){%g,%g,%g};\n",
+ t->getVertex(0)->x(), t->getVertex(0)->y(), t->getVertex(0)->z(),
+ t->getVertex(1)->x(), t->getVertex(1)->y(), t->getVertex(1)->z(),
+ t->getVertex(2)->x(), t->getVertex(2)->y(), t->getVertex(2)->z(),
+ 1., 1., 1.);
+ }
+ fprintf(ftri,"};\n");
+ fclose(ftri);
+ }
+
}
@@ -313,7 +348,31 @@ bool GFaceCompound::trivial() const
}
return false;
}
+//For the conformal map the linear system cannot guarantee there is no folding
+//of triangles
+bool GFaceCompound::checkFolding(std::vector<MVertex*> &ordered) const
+{
+
+ bool has_no_folding = true;
+
+ for(unsigned int i = 0; i < ordered.size()-1; ++i){
+ SPoint3 p1 = coordinates[ordered[i]];
+ SPoint3 p2 = coordinates[ordered[i+1]];
+ int maxSize = (i==0) ? ordered.size()-2: ordered.size()-1;
+ for(unsigned int k = i+2; k < maxSize; ++k){
+ SPoint3 q1 = coordinates[ordered[k]];
+ SPoint3 q2 = coordinates[ordered[k]];
+ double x[2];
+ int inters = intersection_segments (p1,p2,q1,q2,x);
+ if (inters > 0) has_no_folding = false;
+ }
+ }
+
+ if ( !has_no_folding )
+ Msg::Warning("$$$ Folding for compound face %d", this->tag());
+ return has_no_folding;
+}
//check if the discrete harmonic map is correct
//by checking that all the mapped triangles have
@@ -419,52 +478,60 @@ bool GFaceCompound::parametrize() const
if(allNodes.empty()) buildAllNodes();
-// // TEST MULTISCALE LAPLACE ----------------------------------
-// std::vector<MElement*> _elements;
-// for( std::list<GFace*>::const_iterator itt = _compound.begin(); itt != _compound.end(); ++itt)
-// for(unsigned int i = 0; i < (*itt)->triangles.size(); ++i)
-// _elements.push_back((*itt)->triangles[i]);
-// multiscaleLaplace multiLaplace(_elements);
-// printf("ended multiscale \n");
-// exit(1);
- // TEST MULTISCALE LAPLACE ----------------------------------
//Laplace parametrization
//-----------------
if (_mapping == HARMONIC){
- Msg::Debug("Parametrizing surface %d with 'harmonic map'", tag());
+ Msg::Debug("Parametrizing surface %d with 'harmonic map'", tag());
fillNeumannBCS();
parametrize(ITERU,HARMONIC);
parametrize(ITERV,HARMONIC);
}
+ //Multiscale Laplace parametrization
+ //-----------------
+ else if (_mapping == MULTISCALE){
+ std::vector<MElement*> _elements;
+ for( std::list<GFace*>::const_iterator itt = _compound.begin(); itt != _compound.end(); ++itt)
+ for(unsigned int i = 0; i < (*itt)->triangles.size(); ++i)
+ _elements.push_back((*itt)->triangles[i]);
+ multiscaleLaplace multiLaplace(_elements, coordinates);
+ }
//Conformal map parametrization
//-----------------
else if (_mapping == CONFORMAL){
Msg::Debug("Parametrizing surface %d with 'conformal map'", tag());
fillNeumannBCS();
- parametrize_conformal();
+ bool withoutFolding = parametrize_conformal() ;
+ printStuff();
+ if ( withoutFolding == false ){
+ Msg::Warning("$$$ Parametrization switched to harmonic map");
+ parametrize(ITERU,HARMONIC);
+ parametrize(ITERV,HARMONIC);
+ }
}
//Distance function
//-----------------
//compute_distance();
buildOct();
-
+ printStuff();
+
if (!checkOrientation(0)){
- Msg::Info("Parametrization switched to convex combination map");
+ Msg::Info("*** Parametrization switched to convex combination map");
coordinates.clear();
Octree_Delete(oct);
fillNeumannBCS();
parametrize(ITERU,CONVEXCOMBINATION);
parametrize(ITERV,CONVEXCOMBINATION);
+ checkOrientation(0);
buildOct();
}
- printStuff();
+
computeNormals();
- if (!checkAspectRatio()){
+ if (checkAspectRatio() > AR_MAX){
printf("WARNING: geom aspect ratio too high \n");
exit(1);
paramOK = false;
@@ -534,30 +601,28 @@ double GFaceCompound::getSizeH() const
{
SBoundingBox3d bb;
- SOrientedBoundingBox obbox ;
std::vector<SPoint3> vertices;
for(std::set<MVertex *>::iterator itv = allNodes.begin(); itv !=allNodes.end() ; ++itv){
SPoint3 pt((*itv)->x(),(*itv)->y(), (*itv)->z());
vertices.push_back(pt);
bb +=pt;
}
-
- //obbox = SOrientedBoundingBox::buildOBB(vertices);
- //double H = obbox.getMaxSize();
-
double H = norm(SVector3(bb.max(), bb.min()));
+ //SOrientedBoundingBox obbox = SOrientedBoundingBox::buildOBB(vertices);
+ //double H = obbox.getMaxSize();
+
return H;
}
double GFaceCompound::getSizeBB(const std::list<GEdge* > &elist) const
{
- SOrientedBoundingBox obboxD = obb_boundEdges(elist);
- double D = obboxD.getMaxSize();
+ //SOrientedBoundingBox obboxD = obb_boundEdges(elist);
+ //double D = obboxD.getMaxSize();
- //SOrientedBoundingBox bboxD = boundEdges(elist);
- //double D = norm(SVector3(bboxD.max(), bboxD.min()));
+ SBoundingBox3d bboxD = boundEdges(elist);
+ double D = norm(SVector3(bboxD.max(), bboxD.min()));
return D;
@@ -1035,7 +1100,7 @@ void GFaceCompound::parametrize(iterationStep step, typeOfMapping tom) const
}
-void GFaceCompound::parametrize_conformal() const
+bool GFaceCompound::parametrize_conformal() const
{
dofManager<double> myAssembler(_lsys);
@@ -1045,20 +1110,22 @@ void GFaceCompound::parametrize_conformal() const
bool success = orderVertices(_U0, ordered, coords);
if(!success){
Msg::Error("Could not order vertices on boundary");
- return;
+ return false;
}
- MVertex *v1 = ordered[0];
- MVertex *v2 ;
- double maxSize = 0.0;
- for (int i=1; i< ordered.size(); i++){
- MVertex *vi= ordered[i];
- double dist = vi->distance(v1);
- if (dist > maxSize){
- v2 = vi;
- maxSize = dist;
- }
- }
+ MVertex *v1 = ordered[0];
+ MVertex *v2 = ordered[(int)ceil(ordered.size()/2)];
+
+// MVertex *v2 ;
+// double maxSize = 0.0;
+// for (int i=1; i< ordered.size(); i++){
+// MVertex *vi= ordered[i];
+// double dist = vi->distance(v1);
+// if (dist > maxSize){
+// v2 = vi;
+// maxSize = dist;
+// }
+// }
myAssembler.fixVertex(v1, 0, 1, 0);//0
myAssembler.fixVertex(v1, 0, 2, 0);//0
@@ -1127,6 +1194,9 @@ void GFaceCompound::parametrize_conformal() const
_lsys->clear();
+ //check for folding
+ return checkFolding(ordered);
+
}
void GFaceCompound::compute_distance() const
@@ -1611,22 +1681,27 @@ bool GFaceCompound::checkTopology() const
double H = getSizeH();
double D = H;
- if (_interior_loops.size() > 0) D = getSizeBB(_U0);
- int AR = (int) ceil(H/D);
-
+ if (_interior_loops.size() > 0) D = getSizeBB(_U0);
+ int AR = (int) checkAspectRatio();
+ //int AR = (int) ceil(H/D);
+
if (G != 0 || Nb < 1){
correctTopo = false;
- nbSplit = 2; //std::max(G+2, 2);
+ nbSplit = std::max(G+2, 2);
Msg::Info("-----------------------------------------------------------");
Msg::Warning("Wrong topology: Genus=%d, Nb boundaries=%d, AR=%g", G, Nb, H/D);
Msg::Info("*** Split surface %d in %d parts with Multilevel Mesh partitioner", tag(), nbSplit);
}
- else if (G == 0 && AR > 3){
+ else if (G == 0 && AR > AR_MAX){
correctTopo = false;
nbSplit = -2;
Msg::Info("-----------------------------------------------------------");
Msg::Warning("Wrong topology: Genus=%d, Nb boundaries=%d, AR=%d", G, Nb, AR);
Msg::Info("*** Split surface %d in 2 parts with Laplacian Mesh partitioner", tag());
+
+// correctTopo = true;
+// _mapping = MULTISCALE;
+// Msg::Warning("Aspect Ratio (AR=%d) is too high: using multiscale Laplace", AR);
}
else{
Msg::Debug("Correct topology: Genus=%d and Nb boundaries=%d, AR=%g", G, Nb, H/D);
@@ -1636,17 +1711,17 @@ bool GFaceCompound::checkTopology() const
}
-bool GFaceCompound::checkAspectRatio() const
+double GFaceCompound::checkAspectRatio() const
{
//if ((*(_compound.begin()))->geomType() != GEntity::DiscreteSurface)
// return true;
- bool paramOK = true;
if(allNodes.empty()) buildAllNodes();
double limit = 1.e-15;
double areaMin = 1.e15;
+ double area3D = 0.0;
int nb = 0;
std::list<GFace*>::const_iterator it = _compound.begin();
for( ; it != _compound.end() ; ++it){
@@ -1663,6 +1738,7 @@ bool GFaceCompound::checkAspectRatio() const
double p1[3] = {v[1]->x(), v[1]->y(), v[1]->z()};
double p2[3] = {v[2]->x(), v[2]->y(), v[2]->z()};
double a_3D = fabs(triangle_area(p0, p1, p2));
+ area3D += a_3D;
double q0[3] = {it0->second.x(), it0->second.y(), 0.0};
double q1[3] = {it1->second.x(), it1->second.y(), 0.0};
double q2[3] = {it2->second.x(), it2->second.y(), 0.0};
@@ -1672,6 +1748,20 @@ bool GFaceCompound::checkAspectRatio() const
}
}
+ std::list<GEdge*>::const_iterator it0 = _U0.begin();
+ double tot_length = 0;
+ for( ; it0 != _U0.end(); ++it0 ){
+ for(unsigned int i = 0; i < (*it0)->lines.size(); i++ ){
+ MVertex *v0 = (*it0)->lines[i]->getVertex(0);
+ MVertex *v1 = (*it0)->lines[i]->getVertex(1);
+ const double length = sqrt((v0->x() - v1->x()) * (v0->x() - v1->x()) +
+ (v0->y() - v1->y()) * (v0->y() - v1->y()) +
+ (v0->z() - v1->z()) * (v0->z() - v1->z()));
+ tot_length += length;
+ }
+ }
+ double AR = 2*3.14*area3D/(tot_length*tot_length);
+
if (areaMin < limit && nb > 2) {
Msg::Warning("Geometrical aspect ratio too high (a_2D=%g)", areaMin);
SBoundingBox3d bboxH = bounds();
@@ -1679,14 +1769,12 @@ bool GFaceCompound::checkAspectRatio() const
double D = getSizeBB(_U0);
double eta = H/D;
int nbSplit = -2;
- paramOK = true; //false;
}
else {
Msg::Debug("Geometrical aspect ratio is OK :-)");
- paramOK = true;
}
- return paramOK;
+ return AR;
}
@@ -1694,6 +1782,7 @@ void GFaceCompound::coherenceNormals()
{
Msg::Info("Coherence Normals ");
+
std::map<MEdge, std::set<MTriangle*>, Less_Edge > edge2tris;
for(unsigned int i = 0; i < triangles.size(); i++){
MTriangle *t = triangles[i];
diff --git a/Geo/GFaceCompound.h b/Geo/GFaceCompound.h
index 64809a8e3e4b57cd0afbc0367b29e8cfea7b6e50..34928905745a286b7deb120e09f2c5fc09daca20 100644
--- a/Geo/GFaceCompound.h
+++ b/Geo/GFaceCompound.h
@@ -19,6 +19,8 @@
#include "meshGFaceOptimize.h"
#include "linearSystem.h"
+#define AR_MAX 5 //maximal geometrical spacte ratio
+
/*
A GFaceCompound is a model face that is the compound of model faces.
@@ -51,7 +53,7 @@ class Octree;
class GFaceCompound : public GFace {
public:
typedef enum {ITERU=0,ITERV=1,ITERD=2} iterationStep;
- typedef enum {HARMONIC=1,CONFORMAL=2, CONVEXCOMBINATION=3} typeOfMapping;
+ typedef enum {HARMONIC=1,CONFORMAL=2, CONVEXCOMBINATION=3, MULTISCALE=4} typeOfMapping;
typedef enum {UNITCIRCLE, SQUARE} typeOfIsomorphism;
void computeNormals(std::map<MVertex*, SVector3> &normals) const;
protected:
@@ -71,11 +73,12 @@ class GFaceCompound : public GFace {
void buildOct() const ;
void buildAllNodes() const;
void parametrize(iterationStep, typeOfMapping) const;
- void parametrize_conformal() const;
+ bool parametrize_conformal() const;
void compute_distance() const;
bool checkOrientation(int iter) const;
+ bool checkFolding(std::vector<MVertex*> &ordered) const;
void one2OneMap() const;
- bool checkAspectRatio() const;
+ double checkAspectRatio() const;
void computeNormals () const;
void getBoundingEdges();
void getUniqueEdges(std::set<GEdge*> & _unique);
@@ -118,7 +121,7 @@ class GFaceCompound : public GFace {
mutable int nbSplit;
private:
typeOfIsomorphism _type;
- typeOfMapping _mapping;
+ mutable typeOfMapping _mapping;
};
#else
diff --git a/Mesh/meshGFace.cpp b/Mesh/meshGFace.cpp
index e232b67521ba29c7398c46a211c66ea89016ae8d..a88d38b14c75d9f056d8e0c51f03118bfec70605 100644
--- a/Mesh/meshGFace.cpp
+++ b/Mesh/meshGFace.cpp
@@ -1336,8 +1336,6 @@ void partitionAndRemesh(GFaceCompound *gf)
//-----------------------------------------------------
std::list<GFace*> cFaces = gf->getCompounds();
- //PartitionMeshFace(cFaces, options);
-
std::vector<MElement *> elements;
for (std::list<GFace*>::iterator it = cFaces.begin(); it != cFaces.end(); it++)
for(unsigned int j = 0; j < (*it)->getNumMeshElements(); j++)
diff --git a/Mesh/meshPartition.cpp b/Mesh/meshPartition.cpp
index d3ca6b75784a5ba44d84a356b9db25ebefd899a8..07bf7b47880ee7c999b25fefd4500aa94b7f3aeb 100644
--- a/Mesh/meshPartition.cpp
+++ b/Mesh/meshPartition.cpp
@@ -1080,9 +1080,6 @@ void assignPartitionBoundary(GModel *model,
ppv = new partitionVertex(model, -(int)pvertices.size()-1,v2);
pvertices.insert (ppv);
model->add(ppv);
- printf("*** created partitionVertex %d (", ppv->tag());
- for (unsigned int i = 0; i < v2.size(); ++i) printf("%d ", v2[i]);
- printf(")\n");
}
else ppv = *it;
ppv->points.push_back(new MPoint (ve));
@@ -1276,11 +1273,13 @@ void createPartitionFaces(GModel *model, GFaceCompound *gf, int N,
std::list<GFace*>::iterator it = _compound.begin();
for( ; it != _compound.end(); ++it){
+ //printf("tag compound =%d \n", (*it)->tag());
for(unsigned int i = 0; i < (*it)->triangles.size(); ++i){
MTriangle *e = (*it)->triangles[i];
int part = e->getPartition();
for(int j = 0; j < 3; j++){
MVertex *v0 = e->getVertex(j);
+ //printf("v0=%d part=%d\n", v0->getNum(), part); //returns part 0 ???
allNodes[part-1].insert(v0);
}
discreteFaces[part-1]->triangles.push_back(new MTriangle(e->getVertex(0),e->getVertex(1),e->getVertex(2)));
diff --git a/Mesh/multiscalePartition.cpp b/Mesh/multiscalePartition.cpp
index 2df582db91741baeaaab2f05cbaa465a2e4234f7..e437a565eec59c34cf8b59f82d01511549a9b8de 100644
--- a/Mesh/multiscalePartition.cpp
+++ b/Mesh/multiscalePartition.cpp
@@ -134,54 +134,93 @@ static int getGenus (std::vector<MElement *> &elements,
static int getAspectRatio(std::vector<MElement *> &elements,
std::vector<std::vector<MEdge> > &boundaries)
{
- std::set<MVertex*> vs;
- for(unsigned int i = 0; i < elements.size(); i++){
- MElement *e = elements[i];
- for(int j = 0; j < e->getNumVertices(); j++){
- vs.insert(e->getVertex(j));
- }
- }
- SBoundingBox3d bb;
- SOrientedBoundingBox obbox ;
- std::vector<SPoint3> vertices;
- for (std::set<MVertex* >::iterator it = vs.begin(); it != vs.end(); it++){
- SPoint3 pt((*it)->x(),(*it)->y(), (*it)->z());
- vertices.push_back(pt);
- bb += pt;
+
+ double area3D = 0.0;
+ for(unsigned int i = 0; i <elements.size(); ++i){
+ MElement *t = elements[i];
+ std::vector<MVertex *> v(3);
+ for(int k = 0; k < 3; k++) v[k] = t->getVertex(k);
+ double p0[3] = {v[0]->x(), v[0]->y(), v[0]->z()};
+ double p1[3] = {v[1]->x(), v[1]->y(), v[1]->z()};
+ double p2[3] = {v[2]->x(), v[2]->y(), v[2]->z()};
+ double a_3D = fabs(triangle_area(p0, p1, p2));
+ area3D += a_3D;
}
- //obbox = SOrientedBoundingBox::buildOBB(vertices);
- //double H = obbox.getMaxSize();
- double H = norm(SVector3(bb.max(), bb.min()));
-
- double D = H;
- if (boundaries.size() > 0.0 ) D = 0.0;
- for (unsigned int i = 0; i < boundaries.size(); i++){
- std::set<MVertex*> vb;
+ double tot_length = 0.0;
+ for(unsigned int i = 0; i <boundaries.size(); ++i){
std::vector<MEdge> iBound = boundaries[i];
- for (unsigned int j = 0; j < iBound.size(); j++){
- MEdge e = iBound[j];
- vb.insert(e.getVertex(0));
- vb.insert(e.getVertex(1));
- }
- std::vector<SPoint3> vBounds;
- for (std::set<MVertex* >::iterator it = vb.begin(); it != vb.end(); it++){
- SPoint3 pt((*it)->x(),(*it)->y(), (*it)->z());
- vBounds.push_back(pt);
+ double iLength = 0.0;
+ for( unsigned int j = 0; j <iBound.size(); ++j){
+ MVertex *v0 = iBound[j].getVertex(0);
+ MVertex *v1 = iBound[j].getVertex(1);
+ const double length = sqrt((v0->x() - v1->x()) * (v0->x() - v1->x()) +
+ (v0->y() - v1->y()) * (v0->y() - v1->y()) +
+ (v0->z() - v1->z()) * (v0->z() - v1->z()));
+ iLength += length;
}
- SOrientedBoundingBox obboxD = SOrientedBoundingBox::buildOBB(vBounds);
- D = std::max(D, obboxD.getMaxSize());
+ tot_length += iLength;
}
- int AR = (int)ceil(H/D);
+ int AR = 1.0;
+ if (boundaries.size() > 0 ){
+ tot_length /= boundaries.size();
+ AR = (int) ceil(2*3.14*area3D/(tot_length*tot_length));
+ }
+
+// std::set<MVertex*> vs;
+// for(unsigned int i = 0; i < elements.size(); i++){
+// MElement *e = elements[i];
+// for(int j = 0; j < e->getNumVertices(); j++){
+// vs.insert(e->getVertex(j));
+// }
+// }
+// SBoundingBox3d bb;
+// std::vector<SPoint3> vertices;
+// for (std::set<MVertex* >::iterator it = vs.begin(); it != vs.end(); it++){
+// SPoint3 pt((*it)->x(),(*it)->y(), (*it)->z());
+// vertices.push_back(pt);
+// bb += pt;
+// }
+// double H = norm(SVector3(bb.max(), bb.min()));
+
+// //SOrientedBoundingBox obbox = SOrientedBoundingBox::buildOBB(vertices);
+// //double H = obbox.getMaxSize();
+
+// double D = H;
+// if (boundaries.size() > 0 ) D = 0.0;
+// for (unsigned int i = 0; i < boundaries.size(); i++){
+// std::set<MVertex*> vb;
+// std::vector<MEdge> iBound = boundaries[i];
+// for (unsigned int j = 0; j < iBound.size(); j++){
+// MEdge e = iBound[j];
+// vb.insert(e.getVertex(0));
+// vb.insert(e.getVertex(1));
+// }
+// std::vector<SPoint3> vBounds;
+// SBoundingBox3d bb;
+// for (std::set<MVertex* >::iterator it = vb.begin(); it != vb.end(); it++){
+// SPoint3 pt((*it)->x(),(*it)->y(), (*it)->z());
+// vBounds.push_back(pt);
+// bb +=pt;
+// }
+// double iD = norm(SVector3(bb.max(), bb.min()));
+// D = std::max(D, iD);
+
+// //SOrientedBoundingBox obboxD = SOrientedBoundingBox::buildOBB(vBounds);
+// //D = std::max(D, obboxD.getMaxSize());
+// }
+// int AR = (int)ceil(H/D);
return AR;
}
-static void getGenusAndRatio(std::vector<MElement *> &elements, int & genus, int &AR)
+static void getGenusAndRatio(std::vector<MElement *> &elements, int & genus, int &AR, int &NB)
{
std::vector<std::vector<MEdge> > boundaries;
genus = getGenus(elements, boundaries);
+ NB = boundaries.size();
AR = getAspectRatio(elements, boundaries);
+
}
static void partitionRegions(std::vector<MElement*> &elements,
@@ -281,7 +320,8 @@ void multiscalePartition::partition(partitionLevel & level, int nbParts, typeOfP
#if defined(HAVE_METIS) || defined(HAVE_CHACO)
if (method == LAPLACIAN){
- multiscaleLaplace multiLaplace(level.elements, level.recur);
+ std::map<MVertex*, SPoint3> coordinates;
+ multiscaleLaplace multiLaplace(level.elements, coordinates);
}
else if (method == MULTILEVEL){
setNumberOfPartitions(nbParts);
@@ -300,8 +340,8 @@ void multiscalePartition::partition(partitionLevel & level, int nbParts, typeOfP
nextLevel->region = i;
levels.push_back(nextLevel);
- int genus, AR;
- getGenusAndRatio(regions[i], genus, AR);
+ int genus, AR, NB;
+ getGenusAndRatio(regions[i], genus, AR, NB);
printLevel (nextLevel->elements, nextLevel->recur,nextLevel->region);
@@ -316,15 +356,15 @@ void multiscalePartition::partition(partitionLevel & level, int nbParts, typeOfP
nextLevel->recur,nextLevel->region, genus, AR, nbParts);
partition(*nextLevel, nbParts, MULTILEVEL);
}
- else if (genus == 0 && AR > 3 ){
- int nbParts = 2;
- Msg::Info("Mesh partition: level (%d-%d) is ZERO-GENUS (AR=%d) ---> LAPLACIAN partition %d parts",
- nextLevel->recur,nextLevel->region, AR, nbParts);
- partition(*nextLevel, nbParts, LAPLACIAN);
- }
+ else if (genus == 0 && AR > 4 || genus == 0 && NB > 1){
+ int nbParts = 2;
+ Msg::Info("Mesh partition: level (%d-%d) is ZERO-GENUS (AR=%d NB=%d) ---> LAPLACIAN partition %d parts",
+ nextLevel->recur,nextLevel->region, AR, NB, nbParts);
+ partition(*nextLevel, nbParts, LAPLACIAN);
+ }
else {
- Msg::Info("*** Mesh partition: level (%d-%d) is ZERO-GENUS (AR=%d)",
- nextLevel->recur,nextLevel->region, AR);
+ Msg::Info("*** Mesh partition: level (%d-%d) is ZERO-GENUS (AR=%d, NB=%d)",
+ nextLevel->recur,nextLevel->region, AR, NB);
}
}
diff --git a/Solver/multiscaleLaplace.cpp b/Solver/multiscaleLaplace.cpp
index b3ba1f1f7fde39a78d6fc9de98850efd83a43c34..9a53da220f0b71f81397447229252c574b8353a3 100644
--- a/Solver/multiscaleLaplace.cpp
+++ b/Solver/multiscaleLaplace.cpp
@@ -22,6 +22,15 @@
#include "MTriangle.h"
#include "robustPredicates.h"
+//--------------------------------------------------------------
+
+struct sort_pred {
+ bool operator()(const std::pair<SPoint2,multiscaleLaplaceLevel*> &left, const std::pair<SPoint2,multiscaleLaplaceLevel*> &right) {
+ return left.first.x() < right.first.x();
+ }
+};
+
+
//--------------------------------------------------------------
static int intersection_segments_b (SPoint2 &p1, SPoint2 &p2,
SPoint2 &q1, SPoint2 &q2,
@@ -230,12 +239,10 @@ static void recur_compute_centers_ (double R, double a1, double a2,
}
}
- //compute interior_loops
- std::vector<std::vector<MEdge> > boundaries;
- connected_bounds(root->elements, boundaries);
-
//add the center of real holes ...
if (root->children.size()==0 ){//|| boundaries.size()-1 != root->children.size()){
+ std::vector<std::vector<MEdge> > boundaries;
+ connected_bounds(root->elements, boundaries);
for (int i = 0; i < boundaries.size(); i++){
std::vector<MEdge> me = boundaries[i];
SPoint2 c(0.0,0.0);
@@ -251,14 +258,14 @@ static void recur_compute_centers_ (double R, double a1, double a2,
rad = std::max(rad,sqrt ((c.x() - p.x())*(c.x() - p.x())+
(c.y() - p.y())*(c.y() - p.y())));
}
- if (std::abs(rad/root->radius) < 0.8 && abs(rad) < 0.99){
+ if (std::abs(rad/root->radius) < 0.9 && abs(rad) < 0.99){
centers.push_back(std::make_pair(c,zero));
}
}
}
//sort centers
- std::sort(centers.begin(),centers.end());
+ std::sort(centers.begin(),centers.end(), sort_pred());
for (int i=1;i<centers.size()-1;i++){
multiscaleLaplaceLevel* m1 = centers[i-1].second;
@@ -271,7 +278,6 @@ static void recur_compute_centers_ (double R, double a1, double a2,
}
}
-
}
//--------------------------------------------------------------
static void recur_cut_edges_ (multiscaleLaplaceLevel * root,
@@ -492,18 +498,18 @@ static void recur_cut_ (double R, double a1, double a2,
double x[2];
nbIntersect += intersection_segments (centers[j].first,centers[j+1].first,pp,farLeft,x);
}
- if (root->recur != 0){
- if (nbIntersect %2 == 0)
- left.push_back(root->elements[i]);
- else
- right.push_back(root->elements[i]);
- }
- else{
+// if (root->recur != 0){
+// if (nbIntersect %2 == 0)
+// left.push_back(root->elements[i]);
+// else
+// right.push_back(root->elements[i]);
+// }
+// else{
if (nbIntersect %2 != 0)
left.push_back(root->elements[i]);
else
right.push_back(root->elements[i]);
- }
+ // }
}
for (int i=1;i<centers.size()-1;i++){
@@ -511,9 +517,6 @@ static void recur_cut_ (double R, double a1, double a2,
multiscaleLaplaceLevel* m2 = centers[i].second;
multiscaleLaplaceLevel* m3 = centers[i+1].second;
if (m2){
- /*center of the local system is always 0,0
- its relative position to its parent is center
- only 2 angles have to be computed for in and out*/
a1 = myatan2 (centers[i-1].first.y() - m2->center.y() , centers[i-1].first.x() - m2->center.x() );
a2 = myatan2 (centers[i+1].first.y() - m2->center.y() , centers[i+1].first.x() - m2->center.x() );
recur_cut_ (m2->radius, a1, a2, m2, left, right);
@@ -525,6 +528,7 @@ static void recur_cut_ (double R, double a1, double a2,
static void connected_left_right (std::vector<MElement *> &left,
std::vector<MElement *> &right ){
+
//connected left
std::vector<std::vector<MElement*> > subRegionsL;
connectedRegions (left,subRegionsL);
@@ -588,11 +592,9 @@ static void printLevel ( const char* fn,
std::set<MVertex*> vs;
- for (int i=0;i<elements.size();i++){
- for (int j=0;j<elements[i]->getNumVertices();j++){
+ for (int i=0;i<elements.size();i++)
+ for (int j=0;j<elements[i]->getNumVertices();j++)
vs.insert(elements[i]->getVertex(j));
- }
- }
bool binary = false;
FILE *fp = fopen (fn, "w");
@@ -734,7 +736,8 @@ static bool checkOrientation(std::vector<MElement *> &elements,
}
//--------------------------------------------------------------
-multiscaleLaplace::multiscaleLaplace (std::vector<MElement *> &elements, int iPart)
+multiscaleLaplace::multiscaleLaplace (std::vector<MElement *> &elements,
+ std::map<MVertex*, SPoint3> &allCoordinates)
{
//To go through this execute gmsh with the option -optimize_hom
@@ -767,13 +770,19 @@ multiscaleLaplace::multiscaleLaplace (std::vector<MElement *> &elements, int iPa
//Recursively parametrize
root->recur = 0;
root->region = 0;
+ root->scale = 1.0;
parametrize(*root);
+ //fill the coordinates
+ std::vector<double> iScale;
+ std::vector<SPoint2> iCenter;
+ fillCoordinates(*root, allCoordinates, iScale, iCenter);
+
//Compute centers for the cut
recur_compute_centers_ (1.0, M_PI, 0.0, root);
//Partition the mesh in left and right
- cut (elements, iPart);
+ cut (elements);
//---- Testing other cut for partitionning ----
//---- cutEdges and connected_regions ----
@@ -827,6 +836,35 @@ multiscaleLaplace::multiscaleLaplace (std::vector<MElement *> &elements, int iPa
// else left.push_back(_all[i]);
// }
+}
+
+void multiscaleLaplace::fillCoordinates (multiscaleLaplaceLevel & level,
+ std::map<MVertex*, SPoint3> &allCoordinates,
+ std::vector<double> &iScale,
+ std::vector<SPoint2> &iCenter){
+
+ iScale.push_back(level.scale);
+ iCenter.push_back(level.center);
+
+ for(unsigned int i = 0; i < level.elements.size(); ++i){
+ MElement *e = level.elements[i];
+ for(unsigned int j = 0; j<e->getNumVertices(); ++j){
+ MVertex *v = e->getVertex(j);
+ SPoint2 coord = level.coordinates[v];
+ for (int k= iScale.size()-1; k > 0; k--){
+ coord = coord*iScale[k] + iCenter[k];
+ }
+ allCoordinates[v] = SPoint3(coord.x(), coord.y(), 0.0);
+ }
+ }
+
+
+ for (int i=0;i<level.children.size();i++){
+ multiscaleLaplaceLevel* m = level.children[i];
+ fillCoordinates(*m, allCoordinates, iScale, iCenter);
+ }
+
+
}
void multiscaleLaplace::parametrize (multiscaleLaplaceLevel & level){
@@ -864,7 +902,7 @@ void multiscaleLaplace::parametrize (multiscaleLaplaceLevel & level){
MElement *e = level.elements[i];
std::vector<SPoint2> localCoord;
double local_size = localSize(e,solution);
- if (local_size < 1.e-6 * global_size)
+ if (local_size < 1.e-5 * global_size)
tooSmall.push_back(e);
else goodSize.push_back(e);
}
@@ -872,61 +910,62 @@ void multiscaleLaplace::parametrize (multiscaleLaplaceLevel & level){
//Only keep the connected elements vectors goodSize and tooSmall
std::vector<std::vector<MElement*> > regGoodSize;
connectedRegions (goodSize,regGoodSize);
- int index=0;
if (regGoodSize.size() > 0){
+ int index=0;
int maxSize= regGoodSize[0].size();
for (int i=1;i< regGoodSize.size() ; i++){
int size = regGoodSize[i].size();
if(size > maxSize){
- maxSize = size;
index = i;
+ maxSize = size;
}
}
+ goodSize.clear();
+ for (int i=0;i< regGoodSize.size() ; i++){
+ if (i == index) goodSize.insert(goodSize.begin(), regGoodSize[i].begin(), regGoodSize[i].end());
+ else tooSmall.insert(tooSmall.begin(), regGoodSize[i].begin(), regGoodSize[i].end());
+ }
}
- goodSize.clear();
- for (int i=0;i< regGoodSize.size() ; i++){
- if (i == index) goodSize.insert(goodSize.begin(), regGoodSize[i].begin(), regGoodSize[i].end());
- else tooSmall.insert(tooSmall.begin(), regGoodSize[i].begin(), regGoodSize[i].end());
- }
-
+ level.elements.clear();
level.elements = goodSize;
//Add the not too small regions to the level.elements
std::vector<std::vector<MElement*> > regions_, regions ;
+ regions.clear(); regions_.clear();
connectedRegions (tooSmall,regions_);
for (int i=0;i< regions_.size() ; i++){
bool really_small_elements = false;
for (int k=0; k<regions_[i].size() ; k++){
MElement *e = regions_[i][k];
double local_size = localSize(e,solution);
- if (local_size < 1.e-8 * global_size)
+ if (local_size < 1.e-7 * global_size)
really_small_elements = true;
}
- if(really_small_elements && regions_[i].size() > 10) regions.push_back(regions_[i]);
+ if(really_small_elements ) regions.push_back(regions_[i]);
else
level.elements.insert(level.elements.begin(), regions_[i].begin(), regions_[i].end() );
}
- //Fill level.coordinates
+ //Fill level.coordinates
std::set<MVertex*> goodSizev;
for(unsigned int i = 0; i < level.elements.size(); ++i){
MElement *e = level.elements[i];
for(unsigned int j = 0; j<e->getNumVertices(); ++j){
- goodSizev.insert(e->getVertex(j));
- level.coordinates[e->getVertex(j)] = solution[e->getVertex(j)];
+ MVertex *v = e->getVertex(j);
+ goodSizev.insert(v);
+ level.coordinates[v] = solution[v];
}
}
//Save multiscale meshes
char name[245];
- sprintf(name,"multiscale_level%d_region%d_real.msh",level.recur, level.region);
+ sprintf(name,"multiscale_%d_%d_real.msh",level.recur, level.region);
printLevel (name,level.elements,0,2.0);
- sprintf(name,"multiscale_level%d_region%d_param.msh",level.recur, level.region);
+ sprintf(name,"multiscale_%d_%d_param.msh",level.recur, level.region);
printLevel (name,level.elements,&level.coordinates,2.0);
-
//For every small region compute a new parametrization
- Msg::Info("Level (%d-%d): %d connected small regions",level.recur,level.region, regions.size());
+ Msg::Info("Level (%d-%d): %d connected small regions",level.recur, level.region, regions.size());
for (int i=0;i< regions.size() ; i++){
std::set<MVertex*> tooSmallv;
tooSmallv.clear();
@@ -940,7 +979,7 @@ void multiscaleLaplace::parametrize (multiscaleLaplaceLevel & level){
multiscaleLaplaceLevel *nextLevel = new multiscaleLaplaceLevel;
nextLevel->elements = regions[i];
nextLevel->recur = level.recur+1;
- nextLevel->region = i;//Emi add something here
+ nextLevel->region = i;
SBoundingBox3d smallB;
for(std::set<MVertex *>::iterator itv = tooSmallv.begin(); itv !=tooSmallv.end() ; ++itv){
SPoint2 p = solution[*itv];
@@ -1020,21 +1059,23 @@ void multiscaleLaplace::parametrize_method (multiscaleLaplaceLevel & level,
}
}
-void multiscaleLaplace::cut (std::vector<MElement *> &elements, int iPart)
+void multiscaleLaplace::cut (std::vector<MElement *> &elements)
{
std::vector<MElement*> left,right;
recur_cut_ (1.0, M_PI, 0.0, root,left,right);
-
connected_left_right(left, right);
+ if ( elements.size() != left.size()+right.size())
+ Msg::Error("Cutting laplace wrong nb elements (%d) != left + right (%d)", elements.size(), left.size()+right.size());
+
elements.clear();
elements.insert(elements.end(),left.begin(),left.end());
elements.insert(elements.end(),right.begin(),right.end());
- //char name[256];
- //sprintf(name, "laplace_%d.msh", iPart);
- //printLevel (name,elements,0,2.0);
+ printLevel ("multiscale-all.msh",elements, 0,2.0);
+ printLevel ("multiscale-left.msh",left,0,2.0);
+ printLevel ("multiscale-right.msh",right,0,2.0);
}
diff --git a/Solver/multiscaleLaplace.h b/Solver/multiscaleLaplace.h
index 8b456290e2cef5c5732d68f38c5d256e574092cf..0bb5c6eec2fa484c9c2d7799f21593c673193564 100644
--- a/Solver/multiscaleLaplace.h
+++ b/Solver/multiscaleLaplace.h
@@ -5,6 +5,7 @@
#include <map>
#include <set>
#include "SPoint2.h"
+#include "SPoint3.h"
#include "linearSystem.h"
class MElement;
@@ -23,12 +24,17 @@ struct multiscaleLaplaceLevel {
class multiscaleLaplace{
public:
- multiscaleLaplace (std::vector<MElement *> &elements, int iPart=0);
- void cut (std::vector<MElement *> &elements,int iPart=0);
+ multiscaleLaplace (std::vector<MElement *> &elements,
+ std::map<MVertex*, SPoint3> &allCoordinates);
+ void cut (std::vector<MElement *> &elements);
typedef enum {HARMONIC=1,CONFORMAL=2, CONVEXCOMBINATION=3} typeOfMapping;
linearSystem<double> *_lsys;
multiscaleLaplaceLevel* root;
+ void fillCoordinates (multiscaleLaplaceLevel & level,
+ std::map<MVertex*, SPoint3> &allCoordinates,
+ std::vector<double> &iScale,
+ std::vector<SPoint2> &iCenter);
void parametrize (multiscaleLaplaceLevel &);
void parametrize_method (multiscaleLaplaceLevel & level,
std::set<MVertex*> &allNodes,