diff --git a/Geo/GFaceCompound.cpp b/Geo/GFaceCompound.cpp
index 55e516e611643c546655fbcafbd1bc8314d35aa4..278273b0aaa6f73ec7e86bc05e2bee97d3b81d2f 100644
--- a/Geo/GFaceCompound.cpp
+++ b/Geo/GFaceCompound.cpp
@@ -464,10 +464,11 @@ bool GFaceCompound::trivial() const
// For the conformal map the linear system cannot guarantee there is
// no overlapping of triangles
-bool GFaceCompound::checkOverlap() const
+bool GFaceCompound::checkOverlap(std::vector<MVertex *> &vert) const
{
- bool has_no_overlap = true;
-
+ vert.clear();
+ bool has_overlap = false;
+ double EPS = 1.e-2;
for(std::list<std::list<GEdge*> >::const_iterator iloop = _interior_loops.begin();
iloop != _interior_loops.end(); iloop++){
std::list<GEdge*> loop = *iloop;
@@ -488,20 +489,33 @@ bool GFaceCompound::checkOverlap() const
SPoint3 q2 = coordinates[orderedLoop[k+1]];
double x[2];
int inters = intersection_segments (p1,p2,q1,q2,x);
- if (inters > 0){
- has_no_overlap = false;
- break;
+ if (inters && x[1] > EPS && x[1] < 1.-EPS){
+ has_overlap = true;
+ MVertex *v1 = orderedLoop[i];
+ MVertex *v2 = orderedLoop[k];
+ std::set<MVertex *>::iterator it1 = ov.find(v1);
+ std::set<MVertex *>::iterator it2 = ov.find(v2);
+ vert.push_back(v1);
+ vert.push_back(v2);
+ return has_overlap;
+ // if(it1 == ov.end() && it1 == ov.end()){
+ // ov.insert(v1);
+ // ov.insert(v2);
+ // vert.push_back(v1);
+ // vert.push_back(v2);
+ // return has_overlap;
+ // }
}
}
}
}
- if ( !has_no_overlap ) {
+ if (has_overlap ) {
Msg::Debug("Overlap for compound face %d", this->tag());
}
- return has_no_overlap;
+ return has_overlap;
}
@@ -548,7 +562,7 @@ bool GFaceCompound::checkOrientation(int iter) const
}
else if (oriented && iter < iterMax){
Msg::Info("Parametrization is bijective (no flips)");
- printStuff();
+ //printStuff();
}
return oriented;
@@ -622,8 +636,6 @@ bool GFaceCompound::parametrize() const
fillNeumannBCS();
parametrize(ITERU,HARMONIC);
parametrize(ITERV,HARMONIC);
- //parametrize(ITERU,CONVEXCOMBINATION);
- //parametrize(ITERV,CONVEXCOMBINATION);
}
// Multiscale Laplace parametrization
else if (_mapping == MULTISCALE){
@@ -637,13 +649,10 @@ bool GFaceCompound::parametrize() const
else if (_mapping == CONFORMAL){
Msg::Debug("Parametrizing surface %d with 'conformal map'", tag());
fillNeumannBCS();
- bool noOverlap = parametrize_conformal_spectral() ;
- if (!noOverlap){
- Msg::Warning("!!! Overlap: parametrization switched to 'FE conformal' map");
- noOverlap = parametrize_conformal();
- }
- if (!noOverlap || !checkOrientation(0) ){
- Msg::Warning("$$$ Flipping: parametrization switched to 'harmonic' map");
+ bool hasOverlap = parametrize_conformal_spectral();
+ if (hasOverlap || !checkOrientation(0) ){
+ printStuff(33);
+ Msg::Warning("$$$ Overlap or Flipping: parametrization switched to 'harmonic' map");
parametrize(ITERU,HARMONIC);
parametrize(ITERV,HARMONIC);
}
@@ -929,9 +938,6 @@ GFaceCompound::GFaceCompound(GModel *m, int tag, std::list<GFace*> &compound,
MONE = new simpleFunction<double>(-1.0);
for(std::list<GFace*>::iterator it = _compound.begin(); it != _compound.end(); ++it){
- //EMI FIX
- //if ((*it)->tag() == 3) _mapping = CONFORMAL;
-
if(!(*it)){
Msg::Error("Incorrect face in compound surface %d\n", tag);
Msg::Exit(1);
@@ -1139,12 +1145,15 @@ linearSystem<double> *_lsys = 0;
bool GFaceCompound::parametrize_conformal_spectral() const
{
#if !defined(HAVE_PETSC) && !defined(HAVE_SLEPC)
- Msg::Error("Gmsh should be compiled with petsc and slepc for using the conformal map.");
- Msg::Error("Switch to harmonic map or see doc on the wiki for installing petsc and slepc:");
- Msg::Error("https://geuz.org/trac/gmsh/wiki/STLRemeshing (username:gmsh,passwd:gmsh)");
- return false;
- parametrize_conformal();
+
+ //Msg::Error("Gmsh should be compiled with petsc and slepc for using the conformal map.");
+ //Msg::Error("Switch to harmonic map or see doc on the wiki for installing petsc and slepc:");
+ //Msg::Error("https://geuz.org/trac/gmsh/wiki/STLRemeshing (username:gmsh,passwd:gmsh)");
+ Msg::Warning("Slepc not installed: parametrization switched to 'FE conformal' map");
+ return parametrize_conformal(0,NULL,NULL);
+
#else
+
linearSystem <double> *lsysA = new linearSystemPETSc<double>;
linearSystem <double> *lsysB = new linearSystemPETSc<double>;
dofManager<double> myAssembler(lsysA, lsysB);
@@ -1244,18 +1253,28 @@ bool GFaceCompound::parametrize_conformal_spectral() const
coordinates[v] = SPoint3(paramu,paramv,0.0);
k = k+2;
}
-
delete lsysA;
delete lsysB;
+ }
+ else{
+ Msg::Warning("Slepc not converged: parametrization switched to 'FE conformal' map");
+ return parametrize_conformal(0,NULL,NULL);
+ }
- return checkOverlap();
+ std::vector<MVertex *> vert;
+ bool hasOverlap = checkOverlap(vert);
+ if (hasOverlap){
+ Msg::Warning("!!! Overlap: parametrization switched to 'FE conformal' map");
+ printStuff(3);
+ return hasOverlap = parametrize_conformal(0, vert[0], vert[1]);
+ }
+
+ return hasOverlap;
- }
- else return false;
#endif
}
-bool GFaceCompound::parametrize_conformal() const
+bool GFaceCompound::parametrize_conformal(int iter, MVertex *v1, MVertex *v2) const
{
linearSystem<double> *_lsys = 0;
@@ -1273,8 +1292,8 @@ bool GFaceCompound::parametrize_conformal() const
dofManager<double> myAssembler(_lsys);
- MVertex *v1 = _ordered[0];
- MVertex *v2 = _ordered[(int)ceil((double)_ordered.size()/2.)];
+ if (!v1) v1 = _ordered[0];
+ if (!v2) v2 = _ordered[(int)ceil((double)_ordered.size()/2.)];
myAssembler.fixVertex(v1, 0, 1, 1.);
myAssembler.fixVertex(v1, 0, 2, 0.);
myAssembler.fixVertex(v2, 0, 1, -1.);
@@ -1338,8 +1357,18 @@ bool GFaceCompound::parametrize_conformal() const
delete _lsys;
- //check for overlapping triangles
- return checkOverlap();
+
+ //check for overlap and compute new mapping with new pinned vertices
+ std::vector<MVertex *> vert;
+ bool hasOverlap = checkOverlap(vert);
+ if (hasOverlap && iter < 3){
+ printf("**********Loop FE conformal iter (%d) v1=%d v2=%d \n", iter, vert[0]->getNum(), vert[1]->getNum());
+ printStuff(100+iter);
+ return hasOverlap = parametrize_conformal(iter+1, vert[0],vert[1]);
+ }
+ else{
+ return hasOverlap;
+ }
}
diff --git a/Geo/GFaceCompound.h b/Geo/GFaceCompound.h
index 0d37acb896688eacac8a7bdb0d634447f88dbe44..74aef07aad073da7433d332211c2d9df4c72af7a 100644
--- a/Geo/GFaceCompound.h
+++ b/Geo/GFaceCompound.h
@@ -60,6 +60,7 @@ class GFaceCompound : public GFace {
typedef enum {UNITCIRCLE, SQUARE} typeOfIsomorphism;
void computeNormals(std::map<MVertex*, SVector3> &normals) const;
protected:
+ mutable std::set<MVertex *> ov;
mutable GRbf *_rbf;
simpleFunction<double> *ONE;
simpleFunction<double> *MONE;
@@ -83,11 +84,11 @@ class GFaceCompound : public GFace {
void buildOct() const ;
void buildAllNodes() const;
void parametrize(iterationStep, typeOfMapping) const;
- bool parametrize_conformal() const;
+ bool parametrize_conformal(int iter, MVertex *v1, MVertex *v2) const;
bool parametrize_conformal_spectral() const;
void compute_distance() const;
bool checkOrientation(int iter) const;
- bool checkOverlap() const;
+ bool checkOverlap(std::vector<MVertex *> &vert) const;
void one2OneMap() const;
double checkAspectRatio() const;
void computeNormals () const;
@@ -105,7 +106,7 @@ class GFaceCompound : public GFace {
SOrientedBoundingBox obb_boundEdges(const std::list<GEdge* > &elist) const;
void fillNeumannBCS() const;
/* double sumAngles(std::vector<MVertex*> ordered) const; */
-
+
public:
GFaceCompound(GModel *m, int tag, std::list<GFace*> &compound,
std::list<GEdge*> &U0, typeOfMapping typ = HARMONIC, int allowPartition=1);
diff --git a/Geo/GModel.cpp b/Geo/GModel.cpp
index 5c5cc990dffa19ca095fa21500333df34f34546c..52f66e430eda04a5f2eae9594fb1d190e00125ed 100644
--- a/Geo/GModel.cpp
+++ b/Geo/GModel.cpp
@@ -1815,16 +1815,13 @@ GModel *GModel::buildCutGModel(gLevelset *ls, bool cutElem, bool saveTri)
for(int i = 0; i < (int)(sizeof(elements) / sizeof(elements[0])); i++)
cutGM->_storeElementsInEntities(elements[i]);
-
cutGM->_associateEntityWithMeshVertices();
cutGM->_storeVerticesInEntities(vertexMap);
for(int i = 0; i < 4; i++)
cutGM->_storePhysicalTagsInEntities(i, physicals[i]);
- double t2 = Cpu();
-
- Msg::Info("Mesh cutting complete (%g s)", t2 - t1);
+ Msg::Info("Mesh cutting complete (%g s)", Cpu() - t1);
return cutGM;
}
diff --git a/Geo/MElementCut.cpp b/Geo/MElementCut.cpp
index 2aac62c779178c122494dd77dd6b6c0acd0259c9..a438ffba8ae41220d227375428ff34ebf26bfd66 100644
--- a/Geo/MElementCut.cpp
+++ b/Geo/MElementCut.cpp
@@ -853,6 +853,7 @@ static void elementCutMesh(MElement *e, std::vector<const gLevelset *> &RPN,
int lsT = triangles[i]->lsTag();
int c = elements[2].count(lsT) + elements[3].count(lsT) + elements[8].count(lsT);
// the surfaces are cut before the volumes!
+ //bool havePhys = e->
int reg = getBorderTag(lsT, c, newElemTags[2][0], borderElemTags[1]);
int physTag = (!gePhysicals.size()) ? 0 : getBorderTag(lsT, c, newPhysTags[2][0], borderPhysTags[1]);
std::vector<int> phys;
@@ -1209,6 +1210,7 @@ GModel *buildCutMesh(GModel *gm, gLevelset *ls,
return cutGM;
}
+ //CutMesh
newVerticesContainer newVertices;
std::map<int, int> borderElemTags[2]; //map<lsTag,elementary>[line=0,surface=1]
std::map<int, int> borderPhysTags[2]; //map<lstag,physical>[line=0,surface=1]
diff --git a/Geo/gmshLevelset.cpp b/Geo/gmshLevelset.cpp
index c784cf167be50e7af291692cbdcb9a0e14d3aa47..bb30399dc748ce59139de59066682d6e4be8b5c1 100644
--- a/Geo/gmshLevelset.cpp
+++ b/Geo/gmshLevelset.cpp
@@ -2,7 +2,7 @@
#include <queue>
#include <stack>
#include "fullMatrix.h"
-#include "mathEvaluator.h"
+
inline double det3(double d11, double d12, double d13, double d21, double d22, double d23, double d31, double d32, double d33) {
return d11 * (d22 * d33 - d23 * d32) - d21 * (d12 * d33 - d13 * d32) + d31 * (d12 * d23 - d13 * d22);
diff --git a/Geo/gmshLevelset.h b/Geo/gmshLevelset.h
index 9beef8585c70ec2d96d734c9a3b40636c0982878..e95e0091b991ce168b9c1f8470b856552f4c9842 100644
--- a/Geo/gmshLevelset.h
+++ b/Geo/gmshLevelset.h
@@ -41,7 +41,6 @@
#define INTER 14
#define CRACK 15
-
class gLevelset
{
protected:
@@ -121,7 +120,7 @@ class gLevelsetSphere : public gLevelsetPrimitive
protected:
double xc, yc, zc, r;
public:
- gLevelsetSphere (const double &x, const double &y, const double &z, const double &R, int tag) : gLevelsetPrimitive(tag), xc(x), yc(y), zc(z), r(R) {}
+ gLevelsetSphere (const double &x, const double &y, const double &z, const double &R, int tag=1) : gLevelsetPrimitive(tag), xc(x), yc(y), zc(z), r(R) {}
virtual double operator () (const double x, const double y, const double z) const
{return sqrt((xc - x) * (xc - x) + (yc - y) * (yc - y) + (zc - z) * (zc - z)) - r;}
int type() const {return SPHERE;}
@@ -133,11 +132,11 @@ protected:
double a, b, c, d;
public:
// define the plane _a*x+_b*y+_c*z+_d, with outward normal (a,b,c)
- gLevelsetPlane (const double _a, const double _b, const double _c, const double _d, int tag) : gLevelsetPrimitive(tag), a(_a), b(_b), c(_c), d(_d) {}
+ gLevelsetPlane (const double _a, const double _b, const double _c, const double _d, int tag=1) : gLevelsetPrimitive(tag), a(_a), b(_b), c(_c), d(_d) {}
// define the plane passing through the point pt and with outward normal norm
- gLevelsetPlane (const double *pt, const double *norm, int tag);
+ gLevelsetPlane (const double *pt, const double *norm, int tag=1);
// define the plane passing through the 3 points pt1,pt2,pt3 and with outward normal (pt1,pt2)x(pt1,pt3)
- gLevelsetPlane (const double *pt1, const double *pt2, const double *pt3, int tag);
+ gLevelsetPlane (const double *pt1, const double *pt2, const double *pt3, int tag=1);
// copy constructor
gLevelsetPlane(const gLevelsetPlane &lv);
virtual gLevelset * clone() const{return new gLevelsetPlane(*this);}
@@ -174,7 +173,7 @@ protected:
public:
// define the data points
- gLevelsetPoints(fullMatrix<double> &_centers, int tag);
+ gLevelsetPoints(fullMatrix<double> &_centers, int tag=1);
// copy constructor
gLevelsetPoints(const gLevelsetPoints &lv);
virtual gLevelset * clone() const{return new gLevelsetPoints(*this);}
@@ -198,7 +197,7 @@ protected:
public:
gLevelsetQuadric() : gLevelsetPrimitive() {init(); }
- gLevelsetQuadric(int tag) : gLevelsetPrimitive(tag) {init(); }
+ gLevelsetQuadric(int tag=1) : gLevelsetPrimitive(tag) {init(); }
gLevelsetQuadric(const gLevelsetQuadric &);
virtual ~gLevelsetQuadric() {}
double operator () (const double x, const double y, const double z) const;
@@ -208,7 +207,7 @@ public:
class gLevelsetGenCylinder : public gLevelsetQuadric
{
public:
- gLevelsetGenCylinder (const double *pt, const double *dir, const double &R, int tag);
+ gLevelsetGenCylinder (const double *pt, const double *dir, const double &R, int tag=1);
gLevelsetGenCylinder (const gLevelsetGenCylinder& );
virtual gLevelset * clone() const{return new gLevelsetGenCylinder(*this);}
int type() const {return GENCYLINDER;}
@@ -217,7 +216,7 @@ public:
class gLevelsetEllipsoid : public gLevelsetQuadric
{
public:
- gLevelsetEllipsoid (const double *pt, const double *dir, const double &a, const double &b, const double &c, int tag);
+ gLevelsetEllipsoid (const double *pt, const double *dir, const double &a, const double &b, const double &c, int tag=1);
gLevelsetEllipsoid (const gLevelsetEllipsoid& );
virtual gLevelset * clone() const{return new gLevelsetEllipsoid(*this);}
int type() const {return ELLIPS;}
@@ -226,7 +225,7 @@ public:
class gLevelsetCone : public gLevelsetQuadric
{
public:
- gLevelsetCone (const double *pt, const double *dir, const double &angle, int tag);
+ gLevelsetCone (const double *pt, const double *dir, const double &angle, int tag=1);
gLevelsetCone (const gLevelsetCone& );
virtual gLevelset * clone() const{return new gLevelsetCone(*this);}
int type() const {return CONE;}
@@ -236,7 +235,7 @@ class gLevelsetGeneralQuadric : public gLevelsetQuadric
{
public:
gLevelsetGeneralQuadric (const double *pt, const double *dir, const double &x2, const double &y2, const double &z2,
- const double &z, const double &c, int tag);
+ const double &z, const double &c, int tag=1);
gLevelsetGeneralQuadric (const gLevelsetGeneralQuadric& );
virtual gLevelset * clone() const{return new gLevelsetGeneralQuadric(*this);}
int type() const {return QUADRIC;}
@@ -246,7 +245,7 @@ class gLevelsetMathEval: public gLevelsetPrimitive
{
mathEvaluator *_expr;
public:
- gLevelsetMathEval(std::string f, int tag);
+ gLevelsetMathEval(std::string f, int tag=1);
~gLevelsetMathEval(){ if (_expr) delete _expr; }
double operator () (const double x, const double y, const double z) const;
int type() const { return UNKNOWN; }
@@ -259,7 +258,7 @@ class gLevelsetPostView : public gLevelsetPrimitive
int _viewIndex;
OctreePost *_octree;
public:
- gLevelsetPostView(int index, int tag) ;
+ gLevelsetPostView(int index, int tag=1) ;
~gLevelsetPostView(){ if(_octree) delete _octree;}
double operator () (const double x, const double y, const double z) const;
int type() const { return UNKNOWN; }
@@ -417,7 +416,7 @@ public:
// face normal to dir1 and not including pt : tag+4
// face normal to dir1 and including pt : tag+5
gLevelsetBox(const double *pt, const double *dir1, const double *dir2, const double *dir3,
- const double &a, const double &b, const double &c, int tag);
+ const double &a, const double &b, const double &c, int tag=1);
// create a box with the 8 vertices (pt1,...,pt8).
// check if the faces are planar.
// tags of the faces are : face(pt5,pt6,pt7,pt8) : tag+0
@@ -427,7 +426,7 @@ public:
// face(pt2,pt3,pt7,pt6) : tag+4
// face(pt1,pt5,pt8,pt4) : tag+5
gLevelsetBox(const double *pt1, const double *pt2, const double *pt3, const double *pt4,
- const double *pt5, const double *pt6, const double *pt7, const double *pt8, int tag);
+ const double *pt5, const double *pt6, const double *pt7, const double *pt8, int tag=1);
gLevelsetBox(const gLevelsetBox &);
virtual gLevelset * clone() const{return new gLevelsetBox(*this);}
int type() const {return BOX;}
@@ -443,7 +442,7 @@ public:
// tags of the faces are : exterior face : tag+0
// plane face including pt : tag+1
// plane face opposite to pt : tag+2
- gLevelsetCylinder (const double *pt, const double *dir, const double &R, const double &H, int tag);
+ gLevelsetCylinder (const double *pt, const double *dir, const double &R, const double &H, int tag=1);
// create a cylinder : pt is the point in the middle of the cylinder base,
// dir is the direction of the cylinder axis,
// R is the outer radius of the cylinder,
@@ -453,7 +452,7 @@ public:
// plane face including pt : tag+1
// plane face opposite to pt : tag+2
// interior face : tag+3
- gLevelsetCylinder (const double *pt, const double *dir, const double &R, const double &r, const double &H, int tag);
+ gLevelsetCylinder (const double *pt, const double *dir, const double &R, const double &r, const double &H, int tag=1);
gLevelsetCylinder(const gLevelsetCylinder &);
virtual gLevelset * clone() const{return new gLevelsetCylinder(*this);}
int type() const {return CYLINDER;}
@@ -490,7 +489,7 @@ public:
gLevelsetConrod (const double *pt, const double *dir1, const double *dir2,
const double &H1, const double &H2, const double &H3,
const double &R1, const double &r1, const double &R2, const double &r2,
- const double &L1, const double &L2, const double &E, int tag);
+ const double &L1, const double &L2, const double &E, int tag=1);
gLevelsetConrod(const gLevelsetConrod &);
virtual gLevelset * clone() const{return new gLevelsetConrod(*this);}
int type() const {return CONROD;}
diff --git a/Mesh/meshGFace.cpp b/Mesh/meshGFace.cpp
index 11f2eb3caf2c668ef7253e0749102c2dcf3dca6f..45e91986f9dca4a9b3778d08bb244686314a790a 100644
--- a/Mesh/meshGFace.cpp
+++ b/Mesh/meshGFace.cpp
@@ -1910,7 +1910,7 @@ void partitionAndRemesh(GFaceCompound *gf)
int nume = gf->model()->getMaxElementaryNumber(1) + 1;
int numf = gf->model()->getMaxElementaryNumber(2) + 1;
std::vector<discreteFace*> pFaces;
- createPartitionFaces(gf->model(), gf, NF, pFaces);
+ createPartitionFaces(gf->model(), elements, NF, pFaces);
gf->model()->createTopologyFromFaces(pFaces);
@@ -1952,6 +1952,7 @@ void partitionAndRemesh(GFaceCompound *gf)
GFaceCompound *gfc = new GFaceCompound(gf->model(), num_gfc, f_compound, U0,
gf->getTypeOfMapping());
+
gfc->meshAttributes.recombine = gf->meshAttributes.recombine;
gf->model()->add(gfc);
diff --git a/Mesh/meshGFaceOptimize.cpp b/Mesh/meshGFaceOptimize.cpp
index 22b3c2d86ef37d7f76bc2b4055af65e92405fe45..a7f3ef941e81bfbad0bbf76b928919169888c40b 100644
--- a/Mesh/meshGFaceOptimize.cpp
+++ b/Mesh/meshGFaceOptimize.cpp
@@ -1933,7 +1933,7 @@ void recombineIntoQuads(GFace *gf,
void quadsToTriangles(GFace *gf, double minqual)
{
- printf("COUCOU %g\n",minqual);
+
std::vector<MQuadrangle*> qds;
for (int i=0;i<gf->quadrangles.size();i++){
MQuadrangle *q = gf->quadrangles[i];
diff --git a/Mesh/meshPartition.cpp b/Mesh/meshPartition.cpp
index 57717553ddf46e3d686ef3688a77268594e15039..341dcf042af3c11f3cb1c887d651f1c4841f1ada 100644
--- a/Mesh/meshPartition.cpp
+++ b/Mesh/meshPartition.cpp
@@ -1305,7 +1305,7 @@ int CreatePartitionBoundaries(GModel *model, bool createGhostCells)
return 1;
}
-void createPartitionFaces(GModel *model, GFaceCompound *gf, int N,
+void createPartitionFaces(GModel *model, std::vector<MElement *> &elements, int N,
std::vector<discreteFace*> &discreteFaces)
{
#if defined(HAVE_SOLVER)
@@ -1314,29 +1314,20 @@ void createPartitionFaces(GModel *model, GFaceCompound *gf, int N,
std::vector<std::set<MVertex*> > allNodes;
int numMax = model->getMaxElementaryNumber(2) + 1;
for(int i = 0; i < N; i++){
- //printf("*** Created discreteFace %d \n", numMax+i);
discreteFace *face = new discreteFace(model, numMax+i);
discreteFaces.push_back(face);
- model->add(face);//delete this
+ model->add(face); //delete this
std::set<MVertex*> mySet;
allNodes.push_back(mySet);
}
- std::list<GFace*> _compound = gf->getCompounds();
- 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)));
+ for(unsigned int i = 0; i < elements.size(); ++i){
+ MElement *e = elements[i];
+ int part = e->getPartition()-1;
+ for(int j = 0; j < 3; j++){
+ allNodes[part].insert(e->getVertex(j));
}
+ discreteFaces[part]->triangles.push_back(new MTriangle(e->getVertex(0),e->getVertex(1),e->getVertex(2))) ;
}
for(int i = 0; i < N; i++){
@@ -1344,6 +1335,7 @@ void createPartitionFaces(GModel *model, GFaceCompound *gf, int N,
discreteFaces[i]->mesh_vertices.push_back(*it);
}
}
+
#endif
}
diff --git a/Mesh/meshPartition.h b/Mesh/meshPartition.h
index 2f6caf2a2d19977f164fb88313a18cb5325fc42b..1e402ff0cac7c4b6ffd76d80bcef8c5dbc37c890 100644
--- a/Mesh/meshPartition.h
+++ b/Mesh/meshPartition.h
@@ -34,7 +34,7 @@ int CreatePartitionBoundaries(GModel *model, bool createGhostCells);
void splitBoundaryEdges(GModel *model,
std::set<partitionEdge*, Less_partitionEdge> &newEdges);
-void createPartitionFaces(GModel *model, GFaceCompound *gf, int num_parts,
+void createPartitionFaces(GModel *model, std::vector<MElement *> &elements, int num_parts,
std::vector<discreteFace*> &pFaces);
#endif
diff --git a/Mesh/multiscalePartition.cpp b/Mesh/multiscalePartition.cpp
index f777360efc92d3b993ea283f1fe99c6134fdf675..a34d6987c989505dd7cee2214d03e86d06bbce8d 100644
--- a/Mesh/multiscalePartition.cpp
+++ b/Mesh/multiscalePartition.cpp
@@ -234,7 +234,6 @@ static void getGenusAndRatio(std::vector<MElement *> &elements, int & genus, int
AR = getAspectRatio(elements, boundaries);
}
-
static void partitionRegions(std::vector<MElement*> &elements,
std::vector<std::vector<MElement*> > ®ions)
{
@@ -325,7 +324,8 @@ multiscalePartition::multiscalePartition(std::vector<MElement *> &elements,
partition(*level, nbParts, method);
- totalParts = assembleAllPartitions();
+ totalParts = assembleAllPartitions(elements);
+
}
void multiscalePartition::setNumberOfPartitions(int &nbParts)
@@ -365,15 +365,13 @@ void multiscalePartition::partition(partitionLevel & level, int nbParts,
int genus, AR, NB;
getGenusAndRatio(regions[i], genus, AR, NB);
- //printLevel (nextLevel->elements, nextLevel->recur,nextLevel->region);
-
if (genus < 0) {
Msg::Error("Genus partition is negative G=%d!", genus);
exit(1);
}
if (genus != 0 ){
- int nbParts = 2; //std::max(genus+2,2);
+ int nbParts = std::max(genus+2,2);
Msg::Info("Mesh partition: level (%d-%d) is %d-GENUS (AR=%d) ---> MULTILEVEL partition %d parts",
nextLevel->recur,nextLevel->region, genus, AR, nbParts);
partition(*nextLevel, nbParts, MULTILEVEL);
@@ -401,15 +399,17 @@ void multiscalePartition::partition(partitionLevel & level, int nbParts,
#endif
}
-int multiscalePartition::assembleAllPartitions()
+int multiscalePartition::assembleAllPartitions(std::vector<MElement*> & elements)
{
int iPart = 1;
+ elements.clear();
for (unsigned i = 0; i< levels.size(); i++){
partitionLevel *iLevel = levels[i];
if(iLevel->elements.size() > 0){
for (unsigned j = 0; j < iLevel->elements.size(); j++){
MElement *e = iLevel->elements[j];
+ elements.push_back(e);
e->setPartition(iPart);
}
iPart++;
diff --git a/Mesh/multiscalePartition.h b/Mesh/multiscalePartition.h
index e55c24b9d99b084f5da3f1815406ad66866bfcd5..4438cd2bcc444f7ac7d129e9cc12296a8b8720f6 100644
--- a/Mesh/multiscalePartition.h
+++ b/Mesh/multiscalePartition.h
@@ -38,7 +38,7 @@ class multiscalePartition{
public:
multiscalePartition(std::vector<MElement *> &elements, int nbParts,
typeOfPartition method, int allowPartition);
- int assembleAllPartitions();
+ int assembleAllPartitions(std::vector<MElement*> & elements);
void setNumberOfPartitions(int &nbParts);
int getNumberOfParts(){return totalParts;}
};
diff --git a/Parser/Gmsh.tab.cpp b/Parser/Gmsh.tab.cpp
index 7f0a33fdd3d745b225bc91b4eb8f0dbdd86ea5ca..9a4eaec02dcdc2b81701680965c0ed367254166b 100644
--- a/Parser/Gmsh.tab.cpp
+++ b/Parser/Gmsh.tab.cpp
@@ -6371,7 +6371,7 @@ yyreduce:
int t = (int)(yyvsp[(4) - (6)].d);
GModel *GM = GModel::current();
GM->buildCutGModel(FindLevelSet(t)->ls, true, true);
- GM->setVisibility(0);
+ GM->setVisibility(0);
}
else if(!strcmp((yyvsp[(2) - (6)].c), "SplitMesh")){
int t = (int)(yyvsp[(4) - (6)].d);
diff --git a/Solver/eigenSolver.cpp b/Solver/eigenSolver.cpp
index f594831c0b98769699de5414c2021f9fc8c4fcea..4bc489dde202c7e9e5156c658658d409e2ef6dae 100644
--- a/Solver/eigenSolver.cpp
+++ b/Solver/eigenSolver.cpp
@@ -55,7 +55,7 @@ bool eigenSolver::solve(int numEigenValues, std::string which)
// set some default options
_try(EPSSetDimensions(eps, numEigenValues, PETSC_DECIDE, PETSC_DECIDE));
- _try(EPSSetTolerances(eps, 1.e-7, 20));//1.e-6 50
+ _try(EPSSetTolerances(eps, 1.e-7, 20));//1.e-7 20
//_try(EPSSetType(eps, EPSKRYLOVSCHUR)); //default
_try(EPSSetType(eps, EPSARNOLDI));
//_try(EPSSetType(eps, EPSARPACK));
diff --git a/Solver/elasticityTerm.cpp b/Solver/elasticityTerm.cpp
index dc034b2a29458c572071b62e4f0b3615ddc5af69..0654e7e01f94f2b7cadf6699f384dcbae9158f98 100644
--- a/Solver/elasticityTerm.cpp
+++ b/Solver/elasticityTerm.cpp
@@ -79,7 +79,6 @@ void elasticityTerm::elementMatrix(SElement *se, fullMatrix<double> &m) const
}
}
else{
- printf("coucou AAAAAAAAAAAAARGH \n");
se->gradNodalTestFunctions (u, v, w, invjac,GradsT);
for (int j = 0; j < nbSF; j++){
BT(j, 0) = Grads[j][0]; B(0, j) = GradsT[j][0];
diff --git a/Solver/multiscaleLaplace.cpp b/Solver/multiscaleLaplace.cpp
index 97a73f904b65d982ae46b6af4553c0207cd9b954..a7e8d5e03655e7ccb67922070d5fd5ed8da4862b 100644
--- a/Solver/multiscaleLaplace.cpp
+++ b/Solver/multiscaleLaplace.cpp
@@ -31,6 +31,7 @@ struct compareRotatedPoints {
angle = atan2(r.y()-l.y(),r.x()-l.x());
}
bool operator ( ) (const SPoint2 &p1, const SPoint2 &p2) const {
+ //sort from left (x=-1) to right (sin=0), cos(=1)
double x1 = (p1.x()-left.x())*cos(angle) + (p1.y()-left.y())*sin(angle);
double x2 = (p2.x()-left.x())*cos(angle) + (p2.y()-left.y())*sin(angle);
if (x1<x2)return true;
@@ -51,6 +52,38 @@ struct sort_pred {
}
};
+static void sort_centers_dist(std::vector<std::pair<SPoint2,multiscaleLaplaceLevel*> > ¢ers,
+ const SPoint2 &leftP){
+
+ std::vector<std::pair<SPoint2,multiscaleLaplaceLevel*> > myCenters(centers);
+ std::vector<std::pair<SPoint2,multiscaleLaplaceLevel*> > newCenters;
+ SPoint2 lPoint = leftP;
+
+ //printf("size centers =%d \n", myCenters.size());
+ while (!myCenters.empty()){
+ //printf("size centers loop =%d \n", myCenters.size());
+ double dist = 1.e6;
+ int numClosest = 0;
+ for (int i= 0; i < myCenters.size(); i++){
+ SPoint2 c = myCenters[i].first;
+ double distc = sqrt((c.x() - lPoint.x())*(c.x() - lPoint.x())+
+ (c.y() - lPoint.y())*(c.y() - lPoint.y()));
+ //printf("distc =%g \n", distc);
+ if (distc < dist){
+ dist = distc;
+ numClosest = i;
+ }
+ }
+ //printf("numClosest =%d \n", numClosest);
+ lPoint = myCenters[numClosest].first;
+ newCenters.push_back(myCenters[numClosest]);
+ myCenters.erase(myCenters.begin()+numClosest);
+ }
+
+ centers = newCenters;
+
+};
+
//--------------------------------------------------------------
static int intersection_segments_b (SPoint2 &p1, SPoint2 &p2,
@@ -68,6 +101,7 @@ static int intersection_segments_b (SPoint2 &p1, SPoint2 &p2,
double QP2 = robustPredicates::orient2d(Q1,Q2,P2);
}
+
//--------------------------------------------------------------
static void recur_connect (MVertex *v,
std::multimap<MVertex*,MEdge> &v2e,
@@ -245,9 +279,9 @@ static void recur_compute_centers_ (double R, double a1, double a2,
centers.clear();
multiscaleLaplaceLevel* zero = 0;
- SPoint2 PL (R*cos(a1),R*sin(a1));
- SPoint2 PR (R*cos(a2),R*sin(a2));
-
+ double fact = 2.5;
+ SPoint2 PL (fact*R*cos(a1),fact*R*sin(a1));
+ SPoint2 PR (fact*R*cos(a2),fact*R*sin(a2));
std::vector<SPoint2> centersChild;
centersChild.clear();
@@ -286,7 +320,7 @@ static void recur_compute_centers_ (double R, double a1, double a2,
(c.y() - p.y())*(c.y() - p.y())));
}
- //check if the center has not been added
+ //check if the center has not been added
bool newCenter = true;
for (std::vector<SPoint2>::iterator it2 = centersChild.begin(); it2 != centersChild.end(); it2++){
SPoint2 p = *it2;
@@ -307,8 +341,12 @@ static void recur_compute_centers_ (double R, double a1, double a2,
(int)root->children.size());
}
- //sort centers
- std::sort(centers.begin(),centers.end(), sort_pred(PL,PR));
+ //sort centers from left to right
+ //std::sort(centers.begin(),centers.end(), sort_pred(PL,PR));
+
+ //sort from distances
+ sort_centers_dist(centers, PL);
+
centers.insert(centers.begin(), std::make_pair(PL,zero));
centers.push_back(std::make_pair(PR,zero));
@@ -325,25 +363,28 @@ static void recur_compute_centers_ (double R, double a1, double a2,
}
//--------------------------------------------------------------
-static void recur_cut_edges_ (multiscaleLaplaceLevel * root,
- std::multimap<MEdge,MElement*,Less_Edge> &e2e,
+static void recur_cut_edges_ (multiscaleLaplaceLevel *root,
std::map<MEdge,MVertex*,Less_Edge> &cutEdges,
std::set<MVertex*> &cutVertices){
- std::set<MEdge,Less_Edge> allEdges;
- for (std::multimap <MEdge,MElement*,Less_Edge>::iterator it = e2e.begin();
- it != e2e.end() ; ++it){
- allEdges.insert(it->first);
- }
-
const double EPS = 0.001;
+ std::multimap<MEdge,MElement*,Less_Edge> e2e;
+ std::set<MEdge,Less_Edge> allEdges;
+ for (int i=0;i<root->elements.size();++i){
+ for (int j=0;j<root->elements[i]->getNumEdges();j++){
+ e2e.insert(std::make_pair(root->elements[i]->getEdge(j),root->elements[i]));
+ allEdges.insert(root->elements[i]->getEdge(j));
+ }
+ }
+
std::vector<std::pair<SPoint2,multiscaleLaplaceLevel*> > ¢ers = root->cut;
- for (unsigned int i=0;i<centers.size()-1;i++){
+ for (unsigned int i=0;i< centers.size()-1;i++){
SPoint2 p1 = centers[i].first;
SPoint2 p2 = centers[i+1].first;
+ //printf("*************** line p1p2 (%g %g) -- (%g %g) \n",p1.x(),p1.y(),p2.x(),p2.y());
for (std::set <MEdge,Less_Edge>::iterator it = allEdges.begin();
it != allEdges.end() ; ++it){
- if(e2e.count(*it) == 2 && cutEdges.find(*it) == cutEdges.end()){
+ if( cutEdges.find(*it) == cutEdges.end()){//e2e.count(*it) == 2 &&
std::map<MVertex *, SPoint2>::iterator it0 = root->coordinates.find(it->getVertex(0));
std::map<MVertex *, SPoint2>::iterator it1 = root->coordinates.find(it->getVertex(1));
if (it0 != root->coordinates.end() && it1 != root->coordinates.end()){
@@ -352,15 +393,14 @@ static void recur_cut_edges_ (multiscaleLaplaceLevel * root,
double x[2];
int inters = intersection_segments (p1,p2,q1,q2,x);
if (inters && x[1] > EPS && x[1] < 1.-EPS){
- // printf("%g %g -- %g %g -- %g %g -- %g %g\n",p1.x(),p1.y(),p2.x(),p2.y(),q1.x(),q1.y(),q2.x(),q2.y());
- MVertex *newv = new MVertex ((1.-x[1])*it->getVertex(0)->x() + x[1]*it->getVertex(1)->x(),
- (1.-x[1])*it->getVertex(0)->y() + x[1]*it->getVertex(1)->y(),
+ MVertex *newv = new MVertex ((1.-x[1])*it->getVertex(0)->x() + x[1]*it->getVertex(1)->x(),
+ (1.-x[1])*it->getVertex(0)->y() + x[1]*it->getVertex(1)->y(),
(1.-x[1])*it->getVertex(0)->z() + x[1]*it->getVertex(1)->z());
cutEdges[*it] = newv;
root->coordinates[newv] = q1*(1.-x[1]) + q2*x[1] ;
}
- else if (inters && x[1] <= EPS)cutVertices.insert(it->getVertex(0));
- else if (inters && x[1] >= 1.-EPS)cutVertices.insert(it->getVertex(1));
+ else if (inters && x[1] <= EPS) cutVertices.insert(it->getVertex(0));
+ else if (inters && x[1] >= 1.-EPS) cutVertices.insert(it->getVertex(1));
}
}
}
@@ -368,7 +408,7 @@ static void recur_cut_edges_ (multiscaleLaplaceLevel * root,
for (unsigned int i = 0; i < centers.size(); i++){
multiscaleLaplaceLevel* m2 = centers[i].second;
if (m2){
- recur_cut_edges_ (m2,e2e,cutEdges,cutVertices);
+ recur_cut_edges_ (m2,cutEdges,cutVertices);
}
}
}
@@ -384,10 +424,7 @@ static void recur_cut_elements_ (multiscaleLaplaceLevel * root,
for (unsigned int i = 0; i < root->elements.size(); i++){
MVertex *c[3] = {0,0,0};
for (int j=0;j<3;j++){
- MEdge ed = root->elements[i]->getEdge(j);
- // if (cutVertices.find (ed.getVertex(0)) != cutVertices.end() &&
- // cutVertices.find (ed.getVertex(1)) != cutVertices.end() )theCut.insert(ed);
-
+ MEdge ed = root->elements[i]->getEdge(j);
std::map<MEdge,MVertex*,Less_Edge> :: iterator it = cutEdges.find(ed);
if (it != cutEdges.end()){
c[j] = it->second;
@@ -398,58 +435,73 @@ static void recur_cut_elements_ (multiscaleLaplaceLevel * root,
newElements.push_back(new MTriangle (root->elements[i]->getVertex(0),c[0],root->elements[i]->getVertex(2)));
newElements.push_back(new MTriangle (root->elements[i]->getVertex(2),c[0],c[1]));
theCut.insert(MEdge(c[0],c[1]));
- // FIXME should be done !!!!!!
- //delete root->elements[i];
}
else if (c[0] && c[2]){
newElements.push_back(new MTriangle (root->elements[i]->getVertex(0),c[0],c[2]));
- newElements.push_back(new MTriangle (root->elements[i]->getVertex(1),c[0],root->elements[i]->getVertex(2)));
- newElements.push_back(new MTriangle (root->elements[i]->getVertex(2),c[0],c[2]));
+ newElements.push_back(new MTriangle (c[0],root->elements[i]->getVertex(1),root->elements[i]->getVertex(2)));
+ newElements.push_back(new MTriangle (root->elements[i]->getVertex(2),c[2],c[0]));
theCut.insert(MEdge(c[0],c[2]));
- //delete root->elements[i];
}
else if (c[1] && c[2]){
newElements.push_back(new MTriangle (root->elements[i]->getVertex(2),c[2],c[1]));
newElements.push_back(new MTriangle (root->elements[i]->getVertex(0),root->elements[i]->getVertex(1),c[2]));
newElements.push_back(new MTriangle (c[2],root->elements[i]->getVertex(1),c[1]));
theCut.insert(MEdge(c[1],c[2]));
- //delete root->elements[i];
}
else if (c[0]){
newElements.push_back(new MTriangle (root->elements[i]->getVertex(0),c[0],root->elements[i]->getVertex(2)));
newElements.push_back(new MTriangle (root->elements[i]->getVertex(2),c[0],root->elements[i]->getVertex(1)));
- if (cutVertices.find (root->elements[i]->getVertex(0)) != cutVertices.end())
+ if (cutVertices.find (root->elements[i]->getVertex(0)) != cutVertices.end()){
theCut.insert(MEdge(c[0],root->elements[i]->getVertex(0)));
- else if (cutVertices.find (root->elements[i]->getVertex(1)) != cutVertices.end())
+ }
+ else if (cutVertices.find (root->elements[i]->getVertex(1)) != cutVertices.end()) {
theCut.insert(MEdge(c[0],root->elements[i]->getVertex(1)));
- else
+ }
+ else{
theCut.insert(MEdge(c[0],root->elements[i]->getVertex(2)));
- //delete root->elements[i];
+ }
}
else if (c[1]){
newElements.push_back(new MTriangle (root->elements[i]->getVertex(1),c[1],root->elements[i]->getVertex(0)));
newElements.push_back(new MTriangle (root->elements[i]->getVertex(0),c[1],root->elements[i]->getVertex(2)));
- if (cutVertices.find (root->elements[i]->getVertex(1)) != cutVertices.end())
+ if (cutVertices.find (root->elements[i]->getVertex(0)) != cutVertices.end()){
+ theCut.insert(MEdge(c[1],root->elements[i]->getVertex(0)));
+ }
+ else if (cutVertices.find (root->elements[i]->getVertex(1)) != cutVertices.end()) {
theCut.insert(MEdge(c[1],root->elements[i]->getVertex(1)));
- else if (cutVertices.find (root->elements[i]->getVertex(2)) != cutVertices.end())
+ }
+ else{
theCut.insert(MEdge(c[1],root->elements[i]->getVertex(2)));
- else
- theCut.insert(MEdge(c[1],root->elements[i]->getVertex(0)));
- //delete root->elements[i];
+ }
}
else if (c[2]){
- newElements.push_back(new MTriangle (root->elements[i]->getVertex(0),c[2],root->elements[i]->getVertex(1)));
- newElements.push_back(new MTriangle (root->elements[i]->getVertex(1),c[2],root->elements[i]->getVertex(2)));
- if (cutVertices.find (root->elements[i]->getVertex(0)) != cutVertices.end())
+ newElements.push_back(new MTriangle (root->elements[i]->getVertex(0),root->elements[i]->getVertex(1), c[2]));
+ newElements.push_back(new MTriangle (root->elements[i]->getVertex(1),root->elements[i]->getVertex(2), c[2]));
+ if (cutVertices.find (root->elements[i]->getVertex(0)) != cutVertices.end()){
theCut.insert(MEdge(c[2],root->elements[i]->getVertex(0)));
- else if (cutVertices.find (root->elements[i]->getVertex(2)) != cutVertices.end())
- theCut.insert(MEdge(c[2],root->elements[i]->getVertex(2)));
- else
+ }
+ else if (cutVertices.find (root->elements[i]->getVertex(1)) != cutVertices.end()) {
theCut.insert(MEdge(c[2],root->elements[i]->getVertex(1)));
- //delete root->elements[i];
+ }
+ else{
+ theCut.insert(MEdge(c[2],root->elements[i]->getVertex(2)));
+ }
+ }
+ else {
+ newElements.push_back(root->elements[i]);
+ if (cutVertices.find (root->elements[i]->getVertex(0)) != cutVertices.end() &&
+ cutVertices.find (root->elements[i]->getVertex(1)) != cutVertices.end())
+ theCut.insert(MEdge(root->elements[i]->getVertex(0),root->elements[i]->getVertex(1)));
+ if (cutVertices.find (root->elements[i]->getVertex(0)) != cutVertices.end() &&
+ cutVertices.find (root->elements[i]->getVertex(1)) != cutVertices.end())
+ theCut.insert(MEdge(root->elements[i]->getVertex(0),root->elements[i]->getVertex(2)));
+ if (cutVertices.find (root->elements[i]->getVertex(2)) != cutVertices.end() &&
+ cutVertices.find (root->elements[i]->getVertex(1)) != cutVertices.end())
+ theCut.insert(MEdge(root->elements[i]->getVertex(2),root->elements[i]->getVertex(1)));
}
- else newElements.push_back(root->elements[i]);
}
+
+ root->elements.clear();
root->elements = newElements;
_all.insert(_all.end(),newElements.begin(),newElements.end());
for (unsigned int i = 0; i < centers.size(); i++){
@@ -459,22 +511,26 @@ static void recur_cut_elements_ (multiscaleLaplaceLevel * root,
}
}
}
+
//--------------------------------------------------------------
-static void recur_split_ (MElement *e,
- std::multimap<MEdge,MElement*,Less_Edge> &e2e,
- std::set<MElement*> &group,
- std::set<MEdge,Less_Edge> &theCut){
- if (group.find(e) != group.end())return;
- group.insert(e);
+static void recur_leftCut_ (MElement *e,
+ std::multimap<MEdge,MElement*,Less_Edge> &e2e,
+ std::set<MEdge,Less_Edge> &theCut,
+ std::set<MElement*> &leftSet){
+
+ if (leftSet.find(e) != leftSet.end())return;
+ leftSet.insert(e);
+ //printf("insert in left %d \n", e->getNum());
for (int i=0;i<e->getNumEdges();i++){
MEdge ed = e->getEdge(i);
if (theCut.find(ed) == theCut.end()){
for (std::multimap <MEdge,MElement*,Less_Edge>::iterator it = e2e.lower_bound(ed);
it != e2e.upper_bound(ed) ; ++it){
- if (it->second != e)recur_split_ (it->second,e2e,group,theCut);
+ if (it->second != e) recur_leftCut_ (it->second,e2e,theCut, leftSet);
}
}
}
+
}
//--------------------------------------------------------------
@@ -496,7 +552,7 @@ static void recur_connect (const MEdge &e,
}
//--------------------------------------------------------------
-static void connectedRegions (std::vector<MElement*> &elements,
+static void connectedRegions (const std::vector<MElement*> &elements,
std::vector<std::vector<MElement*> > ®ions)
{
std::multimap<MEdge,MElement*,Less_Edge> e2e;
@@ -517,6 +573,28 @@ static void connectedRegions (std::vector<MElement*> &elements,
}
}
//--------------------------------------------------------------
+static void keepConnected(std::vector<MElement*> &goodSize, std::vector<MElement*> &tooSmall){
+
+ std::vector<std::vector<MElement*> > regGoodSize;
+ connectedRegions (goodSize,regGoodSize);
+ if (regGoodSize.size() > 0){
+ int index=0;
+ int maxSize= regGoodSize[0].size();
+ for (unsigned int i=1;i< regGoodSize.size() ; i++){
+ int size = regGoodSize[i].size();
+ if(size > maxSize){
+ index = i;
+ maxSize = size;
+ }
+ }
+ goodSize.clear();
+ for (unsigned 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());
+ }
+ }
+}
+//--------------------------------------------------------------
static void recur_cut_ (double R, double a1, double a2,
multiscaleLaplaceLevel * root,
std::vector<MElement *> &left,
@@ -564,54 +642,8 @@ 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);
- int indexL=0;
-
- if (subRegionsL.size() > 0){
- int maxSize= subRegionsL[0].size();
- for (unsigned int i=1;i< subRegionsL.size() ; i++){
- int size = subRegionsL[i].size();
- if(size > maxSize){
- maxSize = size;
- indexL = i;
- }
- }
- }
-
- left.clear();
- for (unsigned int i=0;i< subRegionsL.size() ; i++){
- if (i == indexL)
- left.insert(left.begin(), subRegionsL[i].begin(), subRegionsL[i].end());
- else
- right.insert(right.begin(), subRegionsL[i].begin(), subRegionsL[i].end());
- }
-
- //connected right
- std::vector<std::vector<MElement*> > subRegionsR;
- connectedRegions (right,subRegionsR);
- int indexR=0;
-
- if (subRegionsR.size() > 0){
- int maxSize= subRegionsR[0].size();
- for (unsigned int i=1;i< subRegionsR.size() ; i++){
- int size = subRegionsR[i].size();
- if(size > maxSize){
- maxSize = size;
- indexR = i;
- }
- }
- }
-
- right.clear();
- for (unsigned int i=0;i< subRegionsR.size() ; i++){
- if (i == indexR)
- right.insert(right.begin(), subRegionsR[i].begin(), subRegionsR[i].end());
- else
- left.insert(left.begin(), subRegionsR[i].begin(), subRegionsR[i].end());
- }
+ keepConnected(left, right);
//assign partitions
for (unsigned int i= 0; i< left.size(); i++)
@@ -621,6 +653,34 @@ static void connected_left_right (std::vector<MElement *> &left,
}
//--------------------------------------------------------------
+static void printCut(std::map<MEdge,MVertex*,Less_Edge> &cutEdges, std::set<MEdge,Less_Edge> &theCut, std::set<MVertex*> cutVertices){
+
+ printf("Writing points.pos \n");
+ std::map<MEdge,MVertex*,Less_Edge>::iterator ite = cutEdges.begin();
+ FILE *f1 = fopen("points.pos","w");
+ fprintf(f1,"View\"\"{\n");
+ for ( ; ite != cutEdges.end();++ite){
+ fprintf(f1,"SP(%g,%g,%g){1.0};\n",ite->second->x(),ite->second->y(),ite->second->z());
+ }
+ std::set<MVertex*>::iterator itv = cutVertices.begin();
+ for ( ; itv != cutVertices.end();++itv){
+ fprintf(f1,"SP(%g,%g,%g){3.0};\n",(*itv)->x(),(*itv)->y(),(*itv)->z());
+ }
+ fprintf(f1,"};\n");
+ fclose(f1);
+
+ printf("Writing edges.pos \n");
+ std::set<MEdge,Less_Edge>::iterator itc = theCut.begin();
+ FILE *f2 = fopen("edges.pos","w");
+ fprintf(f2,"View\"\"{\n");
+ for ( ; itc != theCut.end();++itc){
+ fprintf(f2,"SL(%g,%g,%g,%g,%g,%g){1.0,1.0};\n",itc->getVertex(0)->x(),itc->getVertex(0)->y(),itc->getVertex(0)->z(),
+ itc->getVertex(1)->x(),itc->getVertex(1)->y(),itc->getVertex(1)->z());
+ }
+ fprintf(f2,"};\n");
+ fclose(f2);
+}
+//--------------------------------------------------------------
static void printLevel(const char* fn,
std::vector<MElement *> &elements,
std::map<MVertex*,SPoint2> *coordinates,
@@ -847,62 +907,11 @@ multiscaleLaplace::multiscaleLaplace (std::vector<MElement *> &elements,
//Compute centers for the cut
int nbElems = 0;
recur_compute_centers_ (1.0, M_PI, 0.0, root, nbElems);
- //printf("CENTERS: elements =%d, recur nbElems = %d \n", elements.size(), nbElems);
-
- //Partition the mesh in left and right
- cut(elements);
-
- //---- Testing other cut for partitionning ----
- //---- cutEdges and connected_regions ----
-// std::multimap<MEdge,MElement*,Less_Edge> e2e;
-// for (int i=0;i<elements.size();++i){
-// for (int j=0;j<elements[i]->getNumEdges();j++){
-// e2e.insert(std::make_pair(elements[i]->ge
- //a1 = atan2 (m2->center.y() - centers[i-1].first.y(),m2->center.x() - centers[i-1].first.x());
- //a2 = atan2 (m2->center.y() - centers[i+1].first.y(),m2->center.x() - centers[i+1].first.x()); tEdge(j),elements[i]));
-// }
-// }
-// std::map<MEdge,MVertex*,Less_Edge> cutEdges;
-// std::set<MVertex*> cutVertices;
-// recur_cut_edges_ (root,e2e,cutEdges,cutVertices);
-
-// printf("Writing points.pos \n");
-// std::map<MEdge,MVertex*,Less_Edge>::iterator ite = cutEdges.begin();
-// FILE *f1 = fopen("points.pos","w");
-// fprintf(f1,"View\"\"{\n");
-// for ( ; ite != cutEdges.end();++ite){
-// fprintf(f1,"SP(%g,%g,%g){1.0};\n",ite->second->x(),ite->second->y(),ite->second->z());
-// }
-// fprintf(f1,"};\n");
-// fclose(f1);
-// std::set<MEdge,Less_Edge> theCut;
-// std::vector<MElement*> _all;
-// recur_cut_elements_ (root,cutEdges,cutVertices,theCut,_all);
-
-// printf("Writing edges.pos \n");
-// std::set<MEdge,Less_Edge>::iterator itc = theCut.begin();
-// FILE *f2 = fopen("edges.pos","w");
-// fprintf(f2,"View\"\"{\n");
-// for ( ; itc != theCut.end();++itc){
-// fprintf(f2,"SL(%g,%g,%g,%g,%g,%g){1.0,1.0};\n",itc->getVertex(0)->x(),itc->getVertex(0)->y(),itc->getVertex(0)->z(),
-// itc->getVertex(1)->x(),itc->getVertex(1)->y(),itc->getVertex(1)->z());
-// }
-// fprintf(f2,"};\n");
-// fclose(f2);
-
-// e2e.clear();
-// for (int i=0;i<_all.size();++i){
-// for (int j=0;j<_all[i]->getNumEdges();j++){
-// e2e.insert(std::make_pair(_all[i]->getEdge(j),_all[i]));
-// }
-// }
-// std::set<MElement*> leftSet;
-// recur_split_ (_all[0],e2e,leftSet,theCut);
-// for (int i=0;i<_all.size();i++){
-// if(leftSet.find(_all[i]) == leftSet.end())right.push_back(_all[i]);
-// else left.push_back(_all[i]);
-// }
+ //Split the mesh in left and right
+ //or Cut the mesh in left and right
+ //splitElems(elements);
+ cutElems(elements);
}
@@ -935,8 +944,7 @@ void multiscaleLaplace::fillCoordinates (multiscaleLaplaceLevel & level,
}
-void multiscaleLaplace::parametrize (multiscaleLaplaceLevel & level){
-
+void multiscaleLaplace::parametrize(multiscaleLaplaceLevel & level){
//Compute all nodes for the level
std::set<MVertex*> allNodes;
@@ -970,110 +978,87 @@ 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) //1.e-5
+ if (local_size < 1.e-5*global_size) //1.e-6
tooSmall.push_back(e);
else goodSize.push_back(e);
}
- //Only keep the connected elements vectors goodSize and tooSmall
- std::vector<std::vector<MElement*> > regGoodSize;
- connectedRegions (goodSize,regGoodSize);
- if (regGoodSize.size() > 0){
- int index=0;
- int maxSize= regGoodSize[0].size();
- for (unsigned int i=1;i< regGoodSize.size() ; i++){
- int size = regGoodSize[i].size();
- if(size > maxSize){
- index = i;
- maxSize = size;
- }
- }
- goodSize.clear();
- for (unsigned 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());
- }
- }
+ //Only keep the connected elements vectors goodSize (the rest goes into tooSmall)
+ keepConnected(goodSize, tooSmall);
//Add the not too small regions to the level.elements
std::vector<std::vector<MElement*> > regions_, regions ;
- std::vector<SPoint2> cents;
- std::vector<double> rads;
regions.clear(); regions_.clear();
connectedRegions (tooSmall,regions_);
for (unsigned int i=0;i< regions_.size() ; i++){
bool really_small_elements = false;
- // double totArea = 0.0;
- // SPoint2 center = (0.,0., 0.);
for (unsigned int k=0; k<regions_[i].size() ; k++){
MElement *e = regions_[i][k];
-// SPoint2 p0 = solution[e->getVertex(0)];
-// SPoint2 p1 = solution[e->getVertex(1)];
-// SPoint2 p2 = solution[e->getVertex(2)];
-// SPoint2 ec = (.3*(p0.x()+p1.x()+p2.x()), .3*(p0.y()+p1.y()+p2.y()));
-// center +=ec;
-// double q0[3] = {p0.x(), p0.y(), 0.0};
-// double q1[3] = {p1.x(), p1.y(), 0.0};
-// double q2[3] = {p2.x(), p2.y(), 0.0};
-// double area = fabs(triangle_area(q0, q1, q2));
-// totArea += area;
double local_size = localSize(e,solution);
- if (local_size < 1.e-8 * global_size) //1.e-7
+ if (local_size < 1.e-8*global_size) //1.e-7
really_small_elements = true;
}
- //center *= (1./regions_[i].size());
if(really_small_elements ){
regions.push_back(regions_[i]);
- //cents.push_back(center);
- //rads.push_back(sqrt(totArea/3.14));
}
else
goodSize.insert(goodSize.begin(), regions_[i].begin(), regions_[i].end() );
}
-// //EMI TEST
-// //ensure that small elements are circular patches
-// for (int i=0;i< regions.size() ; i++){
-// SPoint2 c = cents[i];
-// double rad = rads[i];
-// for (std::vector<MElement*>::iterator it = regions[i].begin(); it != regions[i].end(); ++it){
-// MElement *e = *it;
-// SPoint2 p0 = solution[e->getVertex(0)];
-// SPoint2 p1 = solution[e->getVertex(1)];
-// SPoint2 p2 = solution[e->getVertex(2)];
-// SPoint2 ec = (.3*(p0.x()+p1.x()+p2.x()), .3*(p0.y()+p1.y()+p2.y()));
-// double dist = sqrt((ec.x()-c.x())*(ec.x()-c.x()) + (ec.y()-c.y())*(ec.y()-c.y()));
-// std::vector<MElement*>::iterator itp;
-// if (dist > 0.5*rad ) {
-// goodSize.push_back(e);
-// itp = it;
-// it++;
-// regions[i].erase(itp);
-// }
-// }
-// std::vector<std::vector<MElement*> > connRegions ;
-// connectedRegions(regions[i],connRegions);
-// int index=0;
-// int maxSize= connRegions[0].size();
-// for (int j=1;j< connRegions.size() ; j++){
-// int size = connRegions[i].size();
-// if(size > maxSize){
-// maxSize = size;
-// index = j;
-// }
-// }
-// for (int j=0;j< connRegions.size() ; j++){
-// if (j == index){
-// regions[i].clear();
-// regions[i].insert(regions[i].begin(), connRegions[j].begin(), connRegions[j].end());
-// }
-// else{
-// goodSize.insert(goodSize.begin(), connRegions[j].begin(), connRegions[j].end());
-// }
-// }
-// }
- //endTEST EMI
-
+ //check for convex small regions patches
+ for (int i=0;i< regions.size() ; i++){
+ std::vector<MElement*> &elemR = regions[i];
+ v2t_cont adj;
+ buildVertexToElement (elemR,adj);
+ for (std::vector<MElement*>::iterator it = elemR.begin(); it != elemR.end(); ++it){
+ int nbNeigh = 0;
+ MElement *e = *it;
+ v2t_cont :: iterator it0 = adj.find(e->getVertex(0));
+ if(it0 != adj.end()) nbNeigh += it0->second.size();
+ v2t_cont :: iterator it1 = adj.find(e->getVertex(1));
+ if(it1 != adj.end()) nbNeigh += it1->second.size();
+ v2t_cont :: iterator it2 = adj.find(e->getVertex(2));
+ if(it2 != adj.end()) nbNeigh += it2->second.size();
+ std::vector<MElement*>::iterator itp;
+ if (nbNeigh < 12) {
+ goodSize.push_back(e);
+ itp = it;
+ it++;
+ elemR.erase(itp);
+ }
+ }
+ keepConnected(elemR, goodSize);
+ }
+ tooSmall.clear();
+ for (int i=0;i< regions.size() ; i++)
+ tooSmall.insert(tooSmall.begin(), regions[i].begin(), regions[i].end());
+
+ keepConnected(goodSize, tooSmall);
+ regions.clear();
+ connectedRegions (tooSmall,regions);
+
+ //ensure that small elements are circular patches
+ // for (int i=0;i< regions.size() ; i++){
+ // SPoint2 c = cents[i];
+ // double rad = rads[i];
+ // for (std::vector<MElement*>::iterator it = regions[i].begin(); it != regions[i].end(); ++it){
+ // MElement *e = *it;
+ // SPoint2 p0 = solution[e->getVertex(0)];
+ // SPoint2 p1 = solution[e->getVertex(1)];
+ // SPoint2 p2 = solution[e->getVertex(2)];
+ // SPoint2 ec = (.3*(p0.x()+p1.x()+p2.x()), .3*(p0.y()+p1.y()+p2.y()));
+ // double dist = sqrt((ec.x()-c.x())*(ec.x()-c.x()) + (ec.y()-c.y())*(ec.y()-c.y()));
+ // std::vector<MElement*>::iterator itp;
+ // if (dist > 0.5*rad ) {
+ // goodSize.push_back(e);
+ // itp = it;
+ // it++;
+ // regions[i].erase(itp);
+ // }
+ // }
+ // keepConnected(regions[i], goodSize);
+ //}
+
level.elements.clear();
level.elements = goodSize;
@@ -1089,12 +1074,12 @@ void multiscaleLaplace::parametrize (multiscaleLaplaceLevel & level){
}
//Save multiscale meshes
- // std::string name1(level._name+"real.msh");
- // std::string name2(level._name+"param.msh");
- // std::string name3(level._name+"param_small.msh");
- // printLevel (name1.c_str(),level.elements,0,2.2);
- // printLevel (name2.c_str(),level.elements,&level.coordinates,2.2);
- // printLevel_onlysmall (name3.c_str(),level.elements,&level.coordinates,2.2,1.e-15);
+ std::string name1(level._name+"real.msh");
+ std::string name2(level._name+"param.msh");
+ std::string name3(level._name+"param_small.msh");
+ printLevel (name1.c_str(),level.elements,0,2.2);
+ printLevel (name2.c_str(),level.elements,&level.coordinates,2.2);
+ printLevel_onlysmall (name3.c_str(),level.elements,&level.coordinates,2.2,1.e-15);
//For every small region compute a new parametrization
Msg::Info("Level (%d-%d): %d connected small regions",level.recur, level.region, regions.size());
@@ -1208,7 +1193,52 @@ void multiscaleLaplace::parametrize_method (multiscaleLaplaceLevel & level,
delete _lsys;
}
-void multiscaleLaplace::cut(std::vector<MElement *> &elements)
+
+void multiscaleLaplace::cutElems(std::vector<MElement *> &elements)
+{
+
+ std::map<MEdge,MVertex*,Less_Edge> cutEdges;
+ std::set<MEdge,Less_Edge> theCut;
+ std::set<MVertex*> cutVertices;
+ elements.clear();
+
+ recur_cut_edges_ (root, cutEdges,cutVertices);
+ recur_cut_elements_ (root,cutEdges,cutVertices,theCut, elements);
+ printCut(cutEdges, theCut, cutVertices);
+
+ std::multimap<MEdge,MElement*,Less_Edge> e2e;
+ for (int i=0;i<elements.size();++i){
+ for (int j=0;j<elements[i]->getNumEdges();j++){
+ e2e.insert(std::make_pair(elements[i]->getEdge(j),elements[i]));
+ }
+ }
+ std::set<MElement*> leftS;
+ leftS.clear();
+ std::vector<MElement*> left,right;
+ recur_leftCut_ (elements[0], e2e, theCut, leftS);
+
+ for (int i=0;i< elements.size();i++){
+ MElement *e = elements[i];
+ if (leftS.find(e) != leftS.end()) left.push_back(e);
+ else right.push_back(e);
+ }
+
+ connected_left_right(left, right);
+ if (left.size()== 0 || right.size() == 0) {
+ printf("KO size left=%d, right=%d not good (zero elems)\n", left.size(), right.size() );
+ exit(1);
+ }
+
+ elements.clear();
+ elements.insert(elements.end(),left.begin(),left.end());
+ elements.insert(elements.end(),right.begin(),right.end());
+
+ printLevel ("Rootcut-left.msh",left,0,2.2);
+ printLevel ("Rootcut-right.msh",right,0,2.2);
+ printLevel ("Rootcut-all.msh",elements, 0,2.2);
+ //exit(1);
+}
+void multiscaleLaplace::splitElems(std::vector<MElement *> &elements)
{
std::vector<MElement*> left,right;
@@ -1216,22 +1246,22 @@ void multiscaleLaplace::cut(std::vector<MElement *> &elements)
recur_cut_ (1.0, M_PI, 0.0, root,left,right);
connected_left_right(left, right);
- // printLevel ("Rootcut-left.msh",left,0,2.2);
- // printLevel ("Rootcut-right.msh",right,0,2.2);
- // printLevel ("Rootcut-all.msh",elements, 0,2.2);
+ printLevel ("Rootsplit-left.msh",left,0,2.2);
+ printLevel ("Rootsplit-right.msh",right,0,2.2);
+ printLevel ("Rootsplit-all.msh",elements, 0,2.2);
- // printLevel ("Rootcut-left-param.msh",left,&root->coordinates,2.2);
- // printLevel_onlysmall ("Rootcut-left-param10.msh",left,&root->coordinates,2.2,1.e-10);
- // printLevel_onlysmall ("Rootcut-left-param12.msh",left,&root->coordinates,2.2,1.e-12);
- // printLevel_onlysmall ("Rootcut-left-param15.msh",left,&root->coordinates,2.2,1.e-15);
+ printLevel ("Rootsplit-left-param.msh",left,&root->coordinates,2.2);
+ //printLevel_onlysmall ("Rootsplit-left-param10.msh",left,&root->coordinates,2.2,1.e-10);
+ // printLevel_onlysmall ("Rootsplit-left-param12.msh",left,&root->coordinates,2.2,1.e-12);
+ // printLevel_onlysmall ("Rootsplit-left-param15.msh",left,&root->coordinates,2.2,1.e-15);
- // printLevel ("Rootcut-right-param.msh",right,&root->coordinates,2.2);
- // printLevel_onlysmall ("Rootcut-right-param10.msh",right,&root->coordinates,2.2,1.e-10);
- // printLevel_onlysmall ("Rootcut-right-param12.msh",right,&root->coordinates,2.2,1.e-12);
- // printLevel_onlysmall ("Rootcut-right-param15.msh",right,&root->coordinates,2.2,1.e-15);
+ printLevel ("Rootsplit-right-param.msh",right,&root->coordinates,2.2);
+ // printLevel_onlysmall ("Rootsplit-right-param10.msh",right,&root->coordinates,2.2,1.e-10);
+ // printLevel_onlysmall ("Rootsplit-right-param12.msh",right,&root->coordinates,2.2,1.e-12);
+ // printLevel_onlysmall ("Rootsplit-right-param15.msh",right,&root->coordinates,2.2,1.e-15);
- // printLevel_onlysmall ("Rootcut-all-param12.msh",elements,&root->coordinates,2.2,1.e-12);
- // printLevel_onlysmall ("Rootcut-all-param15.msh",elements,&root->coordinates,2.2,1.e-15);
+ // printLevel_onlysmall ("Rootsplit-all-param12.msh",elements,&root->coordinates,2.2,1.e-12);
+ // printLevel_onlysmall ("Rootsplit-all-param15.msh",elements,&root->coordinates,2.2,1.e-15);
if ( elements.size() != left.size()+right.size()) {
Msg::Error("Cutting laplace wrong nb elements (%d) != left + right (%d)", elements.size(), left.size()+right.size());
diff --git a/Solver/multiscaleLaplace.h b/Solver/multiscaleLaplace.h
index 46e0ab77c5f994819a848fcfa6f5cdb5a7e3c4c0..62f57cc8b974d1845d08803bb5fa0a566f25a03f 100644
--- a/Solver/multiscaleLaplace.h
+++ b/Solver/multiscaleLaplace.h
@@ -27,7 +27,8 @@ class multiscaleLaplace{
public:
multiscaleLaplace (std::vector<MElement *> &elements,
std::map<MVertex*, SPoint3> &allCoordinates);
- void cut (std::vector<MElement *> &elements);
+ void cutElems (std::vector<MElement *> &elements);
+ void splitElems (std::vector<MElement *> &elements);
typedef enum {HARMONIC=1,CONFORMAL=2, CONVEXCOMBINATION=3} typeOfMapping;
multiscaleLaplaceLevel* root;
diff --git a/benchmarks/levelset/square.geo b/benchmarks/levelset/square.geo
index 69dfe29f3ca89032c13079fbbd8f3e5f254791b8..e51bf058f167ce1adbde9ae91fbe232dbab8bc55 100644
--- a/benchmarks/levelset/square.geo
+++ b/benchmarks/levelset/square.geo
@@ -27,7 +27,7 @@ Mesh 2;
Levelset Plane (1) = {0,1,0,-7};
//Levelset Point (1) = {{0.1, 2, 0},{3,2,0},{9,2,0},{5,2,0}, {0.1, 2.2, -1},{3,2.5,-1},{9,2,-1},{5,2,-1}, {0.2, 2, 1},{3,2,1},{9,2,1},{5,2,1}};
-Levelset CutMesh {1};
+Levelset CutMeshTri {1};
diff --git a/benchmarks/stl/skull.geo b/benchmarks/stl/skull.geo
index 08a735d84ccb815fde2ecb6eab663970679c3f60..abd757bf369074f59174dd9f85a2aaf0740e4748 100644
--- a/benchmarks/stl/skull.geo
+++ b/benchmarks/stl/skull.geo
@@ -1,5 +1,5 @@
Mesh.RemeshParametrization=0; //(0) harmonic (1) conformal
-Mesh.RemeshAlgorithm=1; //(0) nosplit (1) automatic (2) split metis
+Mesh.RemeshAlgorithm=2; //(0) nosplit (1) automatic (2) split metis
Mesh.CharacteristicLengthFactor=0.1;
Mesh.Algorithm3D = 4; //Frontal (4) Delaunay(1)