diff --git a/Fltk/classificationEditor.cpp b/Fltk/classificationEditor.cpp
index db7ccb344527cfa1cb05774ad8750afeebd6c9a6..d61509385d2b8aab1bcd2f811cbcddb951f974f7 100644
--- a/Fltk/classificationEditor.cpp
+++ b/Fltk/classificationEditor.cpp
@@ -196,29 +196,6 @@ static void class_ok_cb(Fl_Widget *w, void *data)
Msg::StatusBar(3, false, "");
}
-static int maxEdgeNum()
-{
- GModel::eiter it = GModel::current()->firstEdge();
- GModel::eiter ite = GModel::current()->lastEdge();
- int MAXX = 0;
- while(it != ite){
- MAXX = std::max(MAXX, (*it)->tag());
- ++it;
- }
- return MAXX;
-}
-
-static int maxFaceNum()
-{
- GModel::fiter it = GModel::current()->firstFace();
- GModel::fiter ite = GModel::current()->lastFace();
- int MAXX = 0;
- while(it != ite){
- MAXX = std::max(MAXX, (*it)->tag());
- ++it;
- }
- return MAXX;
-}
struct compareMLinePtr
{
@@ -263,7 +240,7 @@ static GEdge *getNewModelEdge(GFace *gf1, GFace *gf2,
std::map<std::pair<int, int>, GEdge*>::iterator it =
newEdges.find(std::make_pair<int, int>(i1, i2));
if(it == newEdges.end()){
- discreteEdge *temporary = new discreteEdge(GModel::current(), maxEdgeNum() + 1, 0, 0);
+ discreteEdge *temporary = new discreteEdge(GModel::current(), GModel::current()->maxEdgeNum() + 1, 0, 0);
printf("add new edge gf1=%d gf2=%d \n", t1, t2);
GModel::current()->add(temporary);
newEdges[std::make_pair<int, int>(i1, i2)] = temporary;
@@ -343,7 +320,7 @@ static void class_color_cb(Fl_Widget* w, void* data)
std::list<MTri3*> ::iterator it = tris.begin();
while(it != tris.end()){
if(!(*it)->isDeleted()){
- discreteFace *temporary = new discreteFace(GModel::current(), maxFaceNum() + 1);
+ discreteFace *temporary = new discreteFace(GModel::current(), GModel::current()->maxFaceNum() + 1);
recurClassify(*it, temporary, lines, reverse);
GModel::current()->add(temporary);
}
@@ -372,7 +349,7 @@ static void class_color_cb(Fl_Widget* w, void* data)
for (std::map<std::pair<int, int>, GEdge*>::iterator it = newEdges.begin() ; it != newEdges.end() ; ++it){
GEdge *ge = it->second;
- printf("NEW edge with tag = %d \n", ge->tag());
+ //printf("NEW edge with tag = %d \n", ge->tag());
std::list<MLine*> segments;
for (unsigned int i = 0; i < ge->lines.size(); i++){
@@ -381,25 +358,20 @@ static void class_color_cb(Fl_Widget* w, void* data)
//for each actual GEdge
while (!segments.empty()) {
-
std::vector<MLine*> myLines;
std::list<MLine*>::iterator it = segments.begin();
-
MVertex *vB = (*it)->getVertex(0);
MVertex *vE = (*it)->getVertex(1);
myLines.push_back(*it);
segments.erase(it);
it++;
-
//printf("***candidate mline %d %d of size %d \n", vB->getNum(), vE->getNum(), segments.size());
for (int i=0; i<2; i++) {
-
for (std::list<MLine*>::iterator it = segments.begin() ; it != segments.end(); ++it){
MVertex *v1 = (*it)->getVertex(0);
MVertex *v2 = (*it)->getVertex(1);
//printf("mline %d %d \n", v1->getNum(), v2->getNum());
-
std::list<MLine*>::iterator itp;
if ( v1 == vE ){
//printf("->push back this mline \n");
@@ -419,39 +391,22 @@ static void class_color_cb(Fl_Widget* w, void* data)
vE = v1;
i=-1;
}
-
if (it == segments.end()) break;
-
}
-
if (vB == vE) break;
-
if (segments.empty()) break;
-
//printf("not found VB=%d vE=%d\n", vB->getNum(), vE->getNum());
MVertex *temp = vB;
vB = vE;
vE = temp;
//printf("not found VB=%d vE=%d\n", vB->getNum(), vE->getNum());
-
}
-
-// printf("************ CANDIDATE NEW EDGE \n");
-// for (std::vector<MLine*>::iterator it = myLines.begin() ; it != myLines.end() ; ++it){
-// MVertex *v1 = (*it)->getVertex(0);
-// MVertex *v2 = (*it)->getVertex(1);
-// printf("Line %d %d \n", v1->getNum(), v2->getNum());
-// }
- GEdge *newGe = new discreteEdge(GModel::current(), maxEdgeNum() + 1, 0, 0);
+ GEdge *newGe = new discreteEdge(GModel::current(), GModel::current()->maxEdgeNum() + 1, 0, 0);
newGe->lines.insert(newGe->lines.end(), myLines.begin(), myLines.end());
GModel::current()->add(newGe);
- printf("create new edge with tag =%d\n", maxEdgeNum());
-
}//end for each actual GEdge
-
}
-
- printf("end new edge with tag \n");
+ //printf("end new edge with tag \n");
for (std::map<std::pair<int, int>, GEdge*>::iterator it = newEdges.begin() ; it != newEdges.end() ; ++it){
GEdge *ge = it->second;
@@ -716,9 +671,9 @@ classificationEditor::classificationEditor()
// saved for the ones that have been saved by the user
// and that will be used for next step
- temporary = new discreteEdge(GModel::current(), maxEdgeNum() + 1, 0, 0);
+ temporary = new discreteEdge(GModel::current(), GModel::current()->maxEdgeNum() + 1, 0, 0);
GModel::current()->add(temporary);
- saved = new discreteEdge(GModel::current(), maxEdgeNum() + 1, 0, 0);
+ saved = new discreteEdge(GModel::current(), GModel::current()->maxEdgeNum() + 1, 0, 0);
GModel::current()->add(saved);
_window->end();
diff --git a/Fltk/classificationEditor.h b/Fltk/classificationEditor.h
index 9c9b32633546ec1b9c4db05ded65a80cc3036fa1..8b51e09ee10bc4896f3efef233eb9a6aa8a533f9 100644
--- a/Fltk/classificationEditor.h
+++ b/Fltk/classificationEditor.h
@@ -10,6 +10,7 @@
#include <FL/Fl_Toggle_Button.H>
#include <FL/Fl_Round_Button.H>
#include <FL/Fl_Window.H>
+#include <FL/Fl_Value_Input.H>
#include "GModel.h"
#include "MElement.h"
#include "ColorTable.h"
diff --git a/Geo/GEdgeCompound.cpp b/Geo/GEdgeCompound.cpp
index 33214081625cc2b07813c049668a63a34290f5d5..74900f2711c358b9eea093be68f022a8d978add8 100644
--- a/Geo/GEdgeCompound.cpp
+++ b/Geo/GEdgeCompound.cpp
@@ -25,7 +25,15 @@ GEdgeCompound::GEdgeCompound(GModel *m, int tag, std::vector<GEdge*> &compound,
for (unsigned int i = 0; i < _compound.size(); i++)
_compound[i]->setCompound(this);
- parametrize();
+ for(std::vector<GEdge*>::iterator it = _compound.begin(); it != _compound.end(); ++it){
+ if(!(*it)){
+ Msg::Error("Incorrect edge in compound edge %d\n", tag);
+ Msg::Exit(1);
+ }
+ }
+
+ parametrize();
+
}
GEdgeCompound::GEdgeCompound(GModel *m, int tag, std::vector<GEdge*> &compound)
diff --git a/Geo/GFace.h b/Geo/GFace.h
index 35404222159af861122b4a1e310068ae2927717e..1eb20a716879e58b39ea2b9663e7aed2466565c7 100644
--- a/Geo/GFace.h
+++ b/Geo/GFace.h
@@ -124,6 +124,7 @@ class GFace : public GEntity
int poincareMesh ();
int genusMesh () { return (poincareMesh() + edgeLoops.size() - 2) / 2; }
virtual int genusGeom ();
+ virtual bool checkTopology() = 0 ;
// return the point on the face corresponding to the given parameter
virtual GPoint point(double par1, double par2) const = 0;
diff --git a/Geo/GFaceCompound.cpp b/Geo/GFaceCompound.cpp
index a7e18bc7449b11f8d0a5a723313c636af680d275..9af2f27800c35e3715e1c28f45cfc7d4835b0a5a 100644
--- a/Geo/GFaceCompound.cpp
+++ b/Geo/GFaceCompound.cpp
@@ -27,6 +27,13 @@
#include "linearSystemCSR.h"
#include "linearSystemFull.h"
#include "linearSystemPETSc.h"
+#include "meshPartitionOptions.h"
+#include "meshPartition.h"
+#include "CreateFile.h"
+#include "Context.h"
+#include "meshGFace.h"
+#include "meshGEdge.h"
+#include "discreteFace.h"
static void fixEdgeToValue(GEdge *ed, double value, dofManager<double, double> &myAssembler)
{
@@ -194,6 +201,10 @@ bool GFaceCompound::trivial() const
bool GFaceCompound::checkOrientation() const
{
+
+ //Only check orientation for stl files (1 patch)
+ if(_compound.size() > 1.0) return true;
+
std::list<GFace*>::const_iterator it = _compound.begin();
double a_old = 0, a_new;
bool oriented = true;
@@ -292,8 +303,8 @@ void GFaceCompound::one2OneMap() const
// printf("Kernel CG: ucg=%g vcg=%g \n", u_cg, v_cg);
// bool testbadCavity = checkCavity(vTri);
-// if(testbadCavity == true ) printf("**** New cavity is KO \n");
-// else printf("-- New cavity is OK \n");
+// if(testbadCavity == true ) printf("**** New cavity is KOÂ \n");
+// else printf("-- New cavity is OKÂ \n");
// for(int i=0; i< vTri.size(); i++){
// MTriangle *t = (MTriangle*) vTri[i];
@@ -327,6 +338,7 @@ void GFaceCompound::parametrize() const
if(trivial()) return;
if(!oct){
+
coordinates.clear();
//Laplace parametrization
@@ -344,10 +356,9 @@ void GFaceCompound::parametrize() const
//-----------------
//compute_distance();
- //checkOrientation();
- computeNormals();
buildOct();
-
+ checkOrientation();
+ computeNormals();
}
}
@@ -357,9 +368,13 @@ void GFaceCompound::getBoundingEdges()
for(std::list<GFace*>::iterator it = _compound.begin(); it != _compound.end(); ++it){
(*it)->setCompound(this);
}
-
- //in case the bounding edges are explicitely given
+
+ std::set<GEdge*> _unique;
+ getUniqueEdges(_unique);
+ std::set<GEdge*>::iterator itf = _unique.begin();
+
if(_U0.size()){
+ //in case the bounding edges are explicitely given
std::list<GEdge*>::const_iterator it = _U0.begin();
for( ; it != _U0.end() ; ++it){
l_edges.push_back(*it);
@@ -370,23 +385,22 @@ void GFaceCompound::getBoundingEdges()
l_edges.push_back(*it);
(*it)->addFace(this);
}
- return;
+ std::list<GEdge*> loop;
+ computeALoop(_unique, loop);
+ while(!_unique.empty()) computeALoop(_unique,loop);
}
-
- std::set<GEdge*> _unique;
- getUniqueEdges(_unique);
- std::set<GEdge*>::iterator itf = _unique.begin();
- for( ; itf != _unique.end(); ++itf){
- l_edges.push_back(*itf);
- (*itf)->addFace(this);
+ else{
+ //in case the bounding edges are NOT explicitely given
+ for( ; itf != _unique.end(); ++itf){
+ l_edges.push_back(*itf);
+ (*itf)->addFace(this);
+ }
+ computeALoop(_unique, _U0);
+ while(!_unique.empty()) computeALoop(_unique,_U1);
}
- computeALoop(_unique, _U0);
- while(!_unique.empty()){
- computeALoop(_unique, _U1);
- printf("%d in unique\n", (int)_unique.size());
- _interior_loops.push_back(_U1);
- }
+ return;
+
}
void GFaceCompound::getUniqueEdges(std::set<GEdge*> &_unique)
@@ -418,10 +432,12 @@ void GFaceCompound::getUniqueEdges(std::set<GEdge*> &_unique)
(*itf)->addFace(this);
}
}
-
void GFaceCompound::computeALoop(std::set<GEdge*> &_unique, std::list<GEdge*> &loop)
{
std::list<GEdge*> _loop;
+
+ if (_unique.empty()) return;
+
while(!_unique.empty()) {
std::set<GEdge*>::iterator it = _unique.begin();
GVertex *vB = (*it)->getBeginVertex();
@@ -472,7 +488,11 @@ void GFaceCompound::computeALoop(std::set<GEdge*> &_unique, std::list<GEdge*> &l
if(found == true) break;
}
+
loop = _loop;
+ _interior_loops.push_back(loop);
+ return;
+
}
GFaceCompound::GFaceCompound(GModel *m, int tag, std::list<GFace*> &compound,
@@ -503,7 +523,7 @@ GFaceCompound::GFaceCompound(GModel *m, int tag, std::list<GFace*> &compound,
Msg::Exit(1);
}
}
-
+
getBoundingEdges();
if(!_U0.size()) _type = UNITCIRCLE;
else if(!_V1.size()) _type = UNITCIRCLE;
@@ -821,7 +841,7 @@ void GFaceCompound::parametrize_conformal() const
}
MVertex *v1 = ordered[0];
- MVertex *v2;
+ MVertex *v2 = ordered[1];
// if (!_interior_loops.empty()){
// std::vector<MVertex*> ordered2;
// bool success2 = orderVertices(_U1, ordered2, coords);
@@ -831,8 +851,6 @@ void GFaceCompound::parametrize_conformal() const
// v2 = ordered[1];
// //v2 = ordered[ordered.size()/2];
// }
-
- v2 = ordered[1];
//printf("Pinned vertex %g %g %g / %g %g %g \n", v1->x(), v1->y(), v1->z(), v2->x(), v2->y(), v2->z());
//exit(1);
@@ -894,7 +912,7 @@ void GFaceCompound::parametrize_conformal() const
coordinates[v] = SPoint3(value1,value2,0.0);
}
- // checkOrientation();
+ checkOrientation();
computeNormals();
buildOct();
_lsys->clear();
@@ -904,10 +922,10 @@ void GFaceCompound::compute_distance() const
{
SBoundingBox3d bbox = GModel::current()->bounds();
double L = norm(SVector3(bbox.max(), bbox.min()));
- double mu = L / 28;
- simpleFunction<double> DIFF(mu * mu), ONE(1.0);
+ double mu = L/28;
+ simpleFunction<double> DIFF(mu * mu), MONE(1.0);
dofManager<double, double> myAssembler(_lsys);
- distanceTerm distance(model(), 1, &DIFF, &ONE);
+ distanceTerm distance(model(), 1, &DIFF, &MONE);
std::vector<MVertex*> ordered;
boundVertices(_U0, ordered);
@@ -1282,32 +1300,145 @@ void GFaceCompound::buildOct() const
nbT = count;
Octree_Arrange(oct);
}
+//Verify topological conditions for computing the harmonic map
+bool GFaceCompound::checkTopology()
+{
+
+ bool correctTopo = true;
+
+ int Nb = _interior_loops.size();
+ int G = genusGeom() ;
+ if( G != 0 || Nb < 1) correctTopo = false;
+
+ if (!correctTopo){
+
+ Msg::Info("--> Wrong topology: Genus=%d and N boundary=%d", G, Nb);
+
+ //Partition the mesh and createTopology fror new faces
+ //-----------------------------------------------------
+ meshPartitionOptions options;
+ int N = std::max(G+1, 2);
+ options = CTX::instance()->partitionOptions;
+ options.num_partitions = N;
+ options.partitioner = 2; //METIS
+ options.algorithm = 1 ;
+ int ier = PartitionMesh(GModel::current(), options);
+ int numv = GModel::current()->maxVertexNum() + 1;
+ int nume = GModel::current()->maxEdgeNum() + 1;
+ int numf = GModel::current()->maxFaceNum() + 1;
+ CreateTopologyFromPartition(GModel::current(), this, N);
+ GModel::current()->createTopologyFromMesh();
+
+ //Remesh new faces (Compound Lines and Compound Surfaces)
+ //-----------------------------------------------------
+
+ discreteFace *gdf = new discreteFace(GModel::current(), numf+2*N);
+ GModel::current()->add(gdf);
+ //printf("*** Remesh face %d in new discreteFace =%d\n", this->tag(), gdf->tag());
+
+ int numEdges = GModel::current()->maxEdgeNum() - nume + 1;
+ //printf("*** Created %d discretEdges \n", numEdges);
+ for (int i=0; i < numEdges; i++){
+ std::vector<GEdge*>e_compound;
+ GEdge *pe = GModel::current()->getEdgeByTag(nume+i);//partition edge
+ int num_gec = GModel::current()->maxEdgeNum() + 1;
+ e_compound.push_back(pe);
+ printf("*** Remeshing discreteEdge %d with CompoundLine %d\n", pe->tag(), num_gec);
+ GEdge *gec = new GEdgeCompound(GModel::current(), num_gec, e_compound);
+
+ meshGEdge mge;
+ mge(gec);//meshing 1D
+
+ GModel::current()->remove(pe);
+ //GModel::current()->add(gec);
+ }
+
+ int numVert = GModel::current()->maxVertexNum() - numv + 1;
+ printf("*** Created %d discreteVert \n", numVert);
+ for (int i=0; i < numEdges; i++){
+ GVertex *pv = GModel::current()->getVertexByTag(numv+i);//partition vertex
+ GModel::current()->remove(pv);
+ }
+
+ std::list<GEdge*> b[4];
+ std::set<MVertex*> all_vertices;
+ for (int i=0; i < N; i++){
+ std::list<GFace*>f_compound;
+ GFace *pf = GModel::current()->getFaceByTag(numf+i);//partition face
+ int num_gfc = numf + N + i ;
+ f_compound.push_back(pf);
+
+ printf("*** Remeshing discreteFace %d with CompoundSurface %d\n", pf->tag(), num_gfc);
+ GFace *gfc = new GFaceCompound(GModel::current(), num_gfc, f_compound,
+ b[0], b[1], b[2], b[3]);
+ meshGFace mgf;
+ mgf(gfc);//meshing 2D
+
+ for (int j=0; j < gfc->triangles.size(); j++){
+ MTriangle *t = gfc->triangles[j];
+ MVertex *v1 = t->getVertex(0);
+ MVertex *v2 = t->getVertex(1);
+ MVertex *v3 = t->getVertex(2);
+ gdf->triangles.push_back(new MTriangle(v1, v2, v3));
+ all_vertices.insert(v1);
+ all_vertices.insert(v2);
+ all_vertices.insert(v3);
+ }
+
+ GModel::current()->remove(pf);
+ //GModel::current()->add(gfc);
+ }
+
+ //Put new mesh in a new discreteFace
+ //-----------------------------------------------------
+ for(std::set<MVertex*>::iterator it = all_vertices.begin(); it != all_vertices.end(); ++it){
+ gdf->mesh_vertices.push_back(*it);
+ }
+
+ printf("*** Mesh of surface %d done by assembly remeshed faces\n", this->tag());
+ gdf->setTag(this->tag());
+ meshStatistics.status = GFace::DONE;
+
+ //CreateOutputFile(CTX::instance()->outputFileName, CTX::instance()->mesh.format);
+ //Msg::Exit(1);
+
+ }
+
+ return correctTopo;
+
+}
-int GFaceCompound::genusGeom()
+int GFaceCompound::genusGeom() const
{
std::list<GFace*>::const_iterator it = _compound.begin();
- std::set<MEdge, Less_Edge> es;
+ std::set<MEdge, Less_Edge> es;
std::set<MVertex*> vs;
int N = 0;
for( ; it != _compound.end() ; ++it){
for(unsigned int i = 0; i < (*it)->triangles.size(); ++i){
N++;
MTriangle *e = (*it)->triangles[i];
- for(int j = 0; j < e->getNumVertices(); j++) vs.insert(e->getVertex(j));
- for(int j = 0; j < e->getNumEdges(); j++) es.insert(e->getEdge(j));
+ for(int j = 0; j < e->getNumVertices(); j++){
+ vs.insert(e->getVertex(j));
+ }
+ for(int j = 0; j < e->getNumEdges(); j++){
+ es.insert(e->getEdge(j));
+ }
}
}
- int poincare = vs.size() - es.size() + N;// = 2g + 2 - b
+ int poincare = vs.size() - es.size() + N;
+
+ return (int)(-poincare + 2 - _interior_loops.size())/2;
- return (poincare - 2 + edgeLoops.size())/2;
}
void GFaceCompound::printStuff() const
{
std::list<GFace*>::const_iterator it = _compound.begin();
- char name1[256], name2[256], name3[256];
+ char name0[256], name1[256], name2[256], name3[256];
char name4[256], name5[256], name6[256];
+ sprintf(name0, "UVAREA-%d.pos", (*it)->tag());
sprintf(name1, "UVX-%d.pos", (*it)->tag());
sprintf(name2, "UVY-%d.pos", (*it)->tag());
sprintf(name3, "UVZ-%d.pos", (*it)->tag());
@@ -1322,6 +1453,7 @@ void GFaceCompound::printStuff() const
// sprintf(name5, "XYZV.pos");
// sprintf(name6, "XYZC.pos");
+ FILE * uva = fopen(name0,"w");
FILE * uvx = fopen(name1,"w");
FILE * uvy = fopen(name2,"w");
FILE * uvz = fopen(name3,"w");
@@ -1329,6 +1461,7 @@ void GFaceCompound::printStuff() const
FILE * xyzv = fopen(name5,"w");
FILE * xyzc = fopen(name6,"w");
+ fprintf(uva,"View \"\"{\n");
fprintf(uvx,"View \"\"{\n");
fprintf(uvy,"View \"\"{\n");
fprintf(uvz,"View \"\"{\n");
@@ -1373,6 +1506,20 @@ void GFaceCompound::printStuff() const
it0->second.z(),it1->second.z(),it2->second.z());
//K1, K2, K3);
+ double p0[3] = {t->getVertex(0)->x(), t->getVertex(0)->y(), t->getVertex(0)->z()};
+ double p1[3] = {t->getVertex(1)->x(), t->getVertex(1)->y(), t->getVertex(1)->z()};
+ double p2[3] = {t->getVertex(2)->x(), t->getVertex(2)->y(), t->getVertex(2)->z()};
+ double a_3D = fabs(triangle_area(p0, p1, p2));
+ 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};
+ double a_2D = fabs(triangle_area(q0, q1, q2));
+ double area = a_3D/a_2D;
+ fprintf(uva,"ST(%g,%g,%g,%g,%g,%g,%g,%g,%g){%g,%g,%g};\n",
+ it0->second.x(), it0->second.y(), 0.0,
+ it1->second.x(), it1->second.y(), 0.0,
+ it2->second.x(), it2->second.y(), 0.0,
+ area, area, area);
fprintf(uvx,"ST(%g,%g,%g,%g,%g,%g,%g,%g,%g){%g,%g,%g};\n",
it0->second.x(), it0->second.y(), 0.0,
it1->second.x(), it1->second.y(), 0.0,
@@ -1390,6 +1537,8 @@ void GFaceCompound::printStuff() const
t->getVertex(0)->z(), t->getVertex(1)->z(), t->getVertex(2)->z());
}
}
+ fprintf(uva,"};\n");
+ fclose(uva);
fprintf(uvx,"};\n");
fclose(uvx);
fprintf(uvy,"};\n");
diff --git a/Geo/GFaceCompound.h b/Geo/GFaceCompound.h
index cde7ea41dd66dfa48986be9a213335a0ab7add2c..8d431236988e73527afa9e672fe42caf722c5c2a 100644
--- a/Geo/GFaceCompound.h
+++ b/Geo/GFaceCompound.h
@@ -95,7 +95,9 @@ class GFaceCompound : public GFace {
void * getNativePtr() const { return 0; }
virtual SPoint2 getCoordinates(MVertex *v) const;
virtual double curvatureMax(const SPoint2 ¶m) const;
- virtual int genusGeom ();
+ virtual int genusGeom () const;
+ virtual bool checkTopology();
+ virtual std::list<GFace*> getCompounds() {return _compound;};
private:
typeOfIsomorphism _type;
typeOfMapping _mapping;
diff --git a/Geo/GModel.cpp b/Geo/GModel.cpp
index da749249b847abedb2f497606fb95743e312af3f..7e89222a25f4d12d1b54683ac428cc61e8e752c7 100644
--- a/Geo/GModel.cpp
+++ b/Geo/GModel.cpp
@@ -128,6 +128,44 @@ bool GModel::empty() const
return vertices.empty() && edges.empty() && faces.empty() && regions.empty();
}
+int GModel::maxVertexNum()
+{
+ viter it = firstVertex();
+ viter itv = lastVertex();
+ int MAXX = 0;
+ while(it != itv){
+ MAXX = std::max(MAXX, (*it)->tag());
+ ++it;
+ }
+ return MAXX;
+
+}
+int GModel::maxEdgeNum()
+{
+ eiter it = firstEdge();
+ eiter ite = lastEdge();
+ int MAXX = 0;
+ while(it != ite){
+ MAXX = std::max(MAXX, (*it)->tag());
+ ++it;
+ }
+ return MAXX;
+
+}
+
+int GModel::maxFaceNum()
+{
+
+ fiter it = firstFace();
+ fiter ite = lastFace();
+ int MAXX = 0;
+ while(it != ite){
+ MAXX = std::max(MAXX, (*it)->tag());
+ ++it;
+ }
+ return MAXX;
+}
+
GRegion *GModel::getRegionByTag(int n) const
{
GEntity tmp((GModel*)this, n);
@@ -982,6 +1020,7 @@ void GModel::createTopologyFromMesh()
std::vector<GEntity*> entities;
getEntities(entities);
+
std::vector<discreteEdge*> discreteEdges;
for(eiter it = firstEdge(); it != lastEdge(); it++)
if((*it)->geomType() == GEntity::DiscreteCurve)
@@ -991,38 +1030,36 @@ void GModel::createTopologyFromMesh()
for(fiter it = firstFace(); it != lastFace(); it++)
if((*it)->geomType() == GEntity::DiscreteSurface)
discreteFaces.push_back((discreteFace*) *it);
+
std::vector<discreteRegion*> discreteRegions;
for(riter it = firstRegion(); it != lastRegion(); it++)
if((*it)->geomType() == GEntity::DiscreteVolume)
discreteRegions.push_back((discreteRegion*) *it);
-
+
Msg::Debug("Creating topology from mesh:");
- Msg::Debug("%d discrete edges", discreteEdges.size());
- Msg::Debug("%d discrete faces", discreteFaces.size());
+ Msg::Debug("%d discrete edges", discreteEdges.size());
+ Msg::Debug("%d discrete faces", discreteFaces.size());
Msg::Debug("%d discrete regions", discreteRegions.size());
// for each discreteRegion, create topology
-
- for (std::vector<discreteRegion*>::iterator it = discreteRegions.begin();
- it != discreteRegions.end(); it++)
+ for (std::vector<discreteRegion*>::iterator it = discreteRegions.begin(); it != discreteRegions.end(); it++)
(*it)->setBoundFaces();
- // for each discreteFace, create topology and if needed create
- // discreteEdges
-
+ // for each discreteFace, create topology and if needed create discreteEdges
int initSizeEdges = discreteEdges.size();
// find boundary edges of each face and put them in a map_edges that
// associates the MEdges with the tags of the adjacent faces
std::map<MEdge, std::vector<int>, Less_Edge > map_edges;
- for (std::vector<discreteFace*>::iterator it = discreteFaces.begin();
- it != discreteFaces.end(); it++)
+ for (std::vector<discreteFace*>::iterator it = discreteFaces.begin(); it != discreteFaces.end(); it++)
(*it)->findEdges(map_edges);
+ //return if no boundary edges (torus, sphere, ..)
+ if (map_edges.empty()) return;
+
// create reverse map, for each face find set of MEdges that are
// candidate for new discrete Edges
- int num = discreteEdges.size() + 1;
std::map<int, std::vector<int> > face2Edges;
while (!map_edges.empty()){
@@ -1089,6 +1126,7 @@ void GModel::createTopologyFromMesh()
}
else{
// for each actual GEdge
+ int num = maxEdgeNum()+1;
while (!myEdges.empty()) {
std::vector<MEdge> myLines;
myLines.clear();
@@ -1130,28 +1168,24 @@ void GModel::createTopologyFromMesh()
discreteEdge *e = new discreteEdge(this, num, 0, 0);
add(e);
discreteEdges.push_back(e);
- std::list<MVertex*> all_vertices;
+
+ std::vector<MVertex*> allV;
for(unsigned int i = 0; i < myLines.size(); i++) {
MVertex *v0 = myLines[i].getVertex(0);
MVertex *v1 = myLines[i].getVertex(1);
e->lines.push_back(new MLine( v0, v1));
- if (std::find(all_vertices.begin(), all_vertices.end(), v0) ==
- all_vertices.end()) all_vertices.push_back(v0);
- if (std::find(all_vertices.begin(), all_vertices.end(), v1) ==
- all_vertices.end()) all_vertices.push_back(v1);
+ if (std::find(allV.begin(), allV.end(), v0) == allV.end()) allV.push_back(v0);
+ if (std::find(allV.begin(), allV.end(), v1) == allV.end()) allV.push_back(v1);
+ v0->setEntity(e);
+ v1->setEntity(e);
}
- e->mesh_vertices.insert(e->mesh_vertices.begin(), all_vertices.begin(),
- all_vertices.end());
-
- for (std::vector<int>::iterator itFace = tagFaces.begin();
- itFace != tagFaces.end(); itFace++) {
- GFace *dFace = getFaceByTag(abs(*itFace));
- for (std::list<MVertex*>::iterator itv = all_vertices.begin();
- itv != all_vertices.end(); itv++) {
- std::vector<MVertex*>::iterator itve =
- std::find(dFace->mesh_vertices.begin(), dFace->mesh_vertices.end(), *itv);
- if (itve != dFace->mesh_vertices.end()) dFace->mesh_vertices.erase(itve);
- (*itv)->setEntity(e);
+ e->mesh_vertices.insert(e->mesh_vertices.begin(), allV.begin(),allV.end());
+
+ for (std::vector<int>::iterator itFace = tagFaces.begin(); itFace != tagFaces.end(); itFace++) {
+ GFace *dFace = getFaceByTag(*itFace);
+ for (std::vector<MVertex*>::iterator itv = allV.begin(); itv != allV.end(); itv++) {
+ std::vector<MVertex*>::iterator itve = std::find(dFace->mesh_vertices.begin(), dFace->mesh_vertices.end(), *itv);
+ if (itve != dFace->mesh_vertices.end()) dFace->mesh_vertices.erase(itve);
}
std::map<int, std::vector<int> >::iterator f2e = face2Edges.find(*itFace);
@@ -1172,19 +1206,20 @@ void GModel::createTopologyFromMesh()
}
// set boundary edges for each face
- for (std::vector<discreteFace*>::iterator it = discreteFaces.begin();
- it != discreteFaces.end(); it++){
+ for (std::vector<discreteFace*>::iterator it = discreteFaces.begin(); it != discreteFaces.end(); it++){
std::map<int, std::vector<int> >::iterator ite = face2Edges.find((*it)->tag());
- std::vector<int> myEdges = ite->second;
- (*it)->setBoundEdges(myEdges);
+ if (ite != face2Edges.end()){
+ std::vector<int> myEdges = ite->second;
+ (*it)->setBoundEdges(myEdges);
+ }
}
// for each discreteEdge, create Topology
- for (std::vector<discreteEdge*>::iterator it = discreteEdges.begin();
- it != discreteEdges.end(); it++){
+ for (std::vector<discreteEdge*>::iterator it = discreteEdges.begin(); it != discreteEdges.end(); it++){
(*it)->createTopo();
(*it)->parametrize();
}
+
}
GModel *GModel::buildCutGModel(gLevelset *ls)
diff --git a/Geo/GModel.h b/Geo/GModel.h
index 96a913fd616e68ae7cc6171b16cba559958b21e0..5812fa461bfc030221880760b0ab5cc06b8d5664 100644
--- a/Geo/GModel.h
+++ b/Geo/GModel.h
@@ -154,6 +154,9 @@ class GModel
int getNumFaces() const { return faces.size(); }
int getNumEdges() const { return edges.size(); }
int getNumVertices() const { return vertices.size(); }
+ int maxVertexNum();
+ int maxEdgeNum();
+ int maxFaceNum();
// quickly check if the model is empty (i.e., if it contains no
// entities)
diff --git a/Geo/OCCFace.h b/Geo/OCCFace.h
index 863d184ea2fe27f6f196e8bc17df096ca1172713..ee170f7f2e898953d7b595307fb3d7b3e4722b0f 100644
--- a/Geo/OCCFace.h
+++ b/Geo/OCCFace.h
@@ -40,6 +40,7 @@ class OCCFace : public GFace {
virtual double curvatures(const SPoint2 ¶m, SVector3 *dirMax, SVector3 *dirMin,
double *curvMax, double *curvMin) const;
surface_params getSurfaceParams() const;
+ virtual bool checkTopology() {return true;};
};
#endif
diff --git a/Geo/discreteEdge.cpp b/Geo/discreteEdge.cpp
index 4caeca19000e2f685b61b105b89452fbc005858e..207b05f30d1721801004e02b502cb3c5cf1825bf 100644
--- a/Geo/discreteEdge.cpp
+++ b/Geo/discreteEdge.cpp
@@ -30,6 +30,7 @@ discreteEdge::discreteEdge(GModel *model, int num, GVertex *_v0, GVertex *_v1)
void discreteEdge::createTopo()
{
+
if(!createdTopo){
orderMLines();
setBoundVertices();
@@ -121,6 +122,17 @@ void discreteEdge::orderMLines()
//lines is now a list of ordered MLines
lines = _m;
+
+ //mesh_vertices
+ mesh_vertices.clear();
+ for (unsigned int i = 0; i < lines.size(); ++i){
+ MVertex *v1 = lines[i]->getVertex(0);
+ MVertex *v2 = lines[i]->getVertex(1);
+ if (std::find(mesh_vertices.begin(), mesh_vertices.end(), v1) == mesh_vertices.end()) mesh_vertices.push_back(v1);
+ if (std::find(mesh_vertices.begin(), mesh_vertices.end(), v2) == mesh_vertices.end()) mesh_vertices.push_back(v2);
+ v1->setEntity(this);
+ v2->setEntity(this);
+ }
//special case reverse orientation
if (lines.size() < 2) return;
@@ -129,6 +141,7 @@ void discreteEdge::orderMLines()
printf("coucou here \n");
for (unsigned int i = 0; i < lines.size(); i++) _orientation[i] = !_orientation[i];
}
+
}
void discreteEdge::setBoundVertices()
@@ -149,21 +162,35 @@ void discreteEdge::setBoundVertices()
}
}
if(!existVertex){
- GVertex *gvB = new discreteVertex(model(),vE->getNum());
+ GVertex *gvB = new discreteVertex(model(), model()->maxVertexNum()+1);
bound_vertices.push_back(gvB);
vE->setEntity(gvB);
gvB->mesh_vertices.push_back(vE);
gvB->points.push_back(new MPoint(gvB->mesh_vertices.back()));
model()->add(gvB);
}
- std::vector<MVertex*>::iterator itve = std::find(mesh_vertices.begin(),
- mesh_vertices.end(), vE);
+ std::vector<MVertex*>::iterator itve = std::find(mesh_vertices.begin(), mesh_vertices.end(), vE);
if (itve != mesh_vertices.end()) mesh_vertices.erase(itve);
+
+ for(GModel::eiter edge = model()->firstEdge(); edge != model()->lastEdge(); edge++){
+ std::vector<MVertex*>::iterator itve = std::find((*edge)->mesh_vertices.begin(), (*edge)->mesh_vertices.end(), vE);
+ if (itve != (*edge)->mesh_vertices.end()) (*edge)->mesh_vertices.erase(itve);
+ }
+ for(GModel::fiter face = model()->firstFace(); face != model()->lastFace(); face++){
+ std::vector<MVertex*>::iterator itve = std::find((*face)->mesh_vertices.begin(), (*face)->mesh_vertices.end(), vE);
+ if (itve != (*face)->mesh_vertices.end()) (*face)->mesh_vertices.erase(itve);
+ }
+ for(GModel::riter reg = model()->firstRegion(); reg != model()->lastRegion(); reg++){
+ std::vector<MVertex*>::iterator itve = std::find((*reg)->mesh_vertices.begin(), (*reg)->mesh_vertices.end(), vE);
+ if (itve != (*reg)->mesh_vertices.end()) (*reg)->mesh_vertices.erase(itve);
+ }
+
}
v0 = bound_vertices[0];
v1 = bound_vertices[1];
}
else if (boundv.size() == 0){
+ //printf("closed loop for discrete Edge =%d \n", this->tag());
GVertex* bound_vertex;
std::vector<MLine*>::const_iterator it = lines.begin();
MVertex* vE = (*it)->getVertex(0);
@@ -172,25 +199,41 @@ void discreteEdge::setBoundVertices()
if ((*point)->mesh_vertices[0]->getNum() == vE->getNum()){
bound_vertex = (*point);
existVertex = true;
+ //printf("vertex exist (%g,%g,%g)\n", vE->x(), vE->y(), vE->z());
break;
}
}
if(!existVertex){
- GVertex *gvB = new discreteVertex(model(),vE->getNum());
+ GVertex *gvB = new discreteVertex(model(), model()->maxVertexNum()+1);
bound_vertex = gvB;
+ vE->setEntity(gvB);
gvB->mesh_vertices.push_back(vE);
gvB->points.push_back(new MPoint(gvB->mesh_vertices.back()));
model()->add(gvB);
}
- std::vector<MVertex*>::iterator itve = std::find(mesh_vertices.begin(),
- mesh_vertices.end(), vE);
- if (itve != mesh_vertices.end()) mesh_vertices.erase(itve);
+ std::vector<MVertex*>::iterator itve = std::find(mesh_vertices.begin(),mesh_vertices.end(), vE);
+ if (itve != mesh_vertices.end()) mesh_vertices.erase(itve);
+
+ for(GModel::eiter edge = model()->firstEdge(); edge != model()->lastEdge(); edge++){
+ std::vector<MVertex*>::iterator itve = std::find((*edge)->mesh_vertices.begin(), (*edge)->mesh_vertices.end(), vE);
+ if (itve != (*edge)->mesh_vertices.end()) (*edge)->mesh_vertices.erase(itve);
+ }
+ for(GModel::fiter face = model()->firstFace(); face != model()->lastFace(); face++){
+ std::vector<MVertex*>::iterator itve = std::find((*face)->mesh_vertices.begin(), (*face)->mesh_vertices.end(), vE);
+ if (itve != (*face)->mesh_vertices.end()) (*face)->mesh_vertices.erase(itve);
+ }
+ for(GModel::riter reg = model()->firstRegion(); reg != model()->lastRegion(); reg++){
+ std::vector<MVertex*>::iterator itve = std::find((*reg)->mesh_vertices.begin(), (*reg)->mesh_vertices.end(), vE);
+ if (itve != (*reg)->mesh_vertices.end()) (*reg)->mesh_vertices.erase(itve);
+ }
+
v0 = bound_vertex;
v1 = bound_vertex;
}
v0->addEdge(this);
v1->addEdge(this);
+
}
/*
@@ -218,8 +261,7 @@ void discreteEdge::parametrize()
if (_orientation[i] == 1 ) vR = lines[i]->getVertex(1);
else vR = lines[i]->getVertex(0);
int param = i+1;
- MVertex *vNEW = new MEdgeVertex(vR->x(),vR->y(),vR->z(), this, param,
- -1., vR->getNum());
+ MVertex *vNEW = new MEdgeVertex(vR->x(),vR->y(),vR->z(), this, param, -1., vR->getNum());
old2new.insert(std::make_pair(vR,vNEW));
newVertices.push_back(vNEW);
newLines.push_back(new MLine(vL, vNEW));
@@ -231,6 +273,7 @@ void discreteEdge::parametrize()
_orientation[i] = 1;
mesh_vertices = newVertices;
+
deleteVertexArrays();
lines.clear();
lines = newLines;
diff --git a/Geo/discreteFace.cpp b/Geo/discreteFace.cpp
index ac9995b337b452484d42a677e1b4d409cdc1e292..1c58800b397c32a8c87cc4ebddf1687f8283872a 100644
--- a/Geo/discreteFace.cpp
+++ b/Geo/discreteFace.cpp
@@ -55,9 +55,8 @@ void discreteFace::findEdges(std::map<MEdge, std::vector<int>, Less_Edge> &map_e
void discreteFace::setBoundEdges(std::vector<int> tagEdges)
{
- for (std::vector<int>::iterator it = tagEdges.begin();
- it != tagEdges.end(); it++){
- GEdge *ge = GModel::current()->getEdgeByTag(abs(*it));
+ for (std::vector<int>::iterator it = tagEdges.begin();it != tagEdges.end(); it++){
+ GEdge *ge = GModel::current()->getEdgeByTag(*it);
l_edges.push_back(ge);
l_dirs.push_back(1);
ge->addFace(this);
diff --git a/Geo/discreteFace.h b/Geo/discreteFace.h
index 87ed1db371b2153eefba5257e1fe3318e3e3baab..f210ebf2a3a83c62ee9d48fd7eec68f60b039105 100644
--- a/Geo/discreteFace.h
+++ b/Geo/discreteFace.h
@@ -24,6 +24,7 @@ class discreteFace : public GFace {
SVector3 *dudu, SVector3 *dvdv, SVector3 *dudv) const;
void setBoundEdges(std::vector<int> tagEdges);
void findEdges(std::map<MEdge, std::vector<int>, Less_Edge > &map_edges);
+ virtual bool checkTopology() {return true;};
};
#endif
diff --git a/Geo/gmshFace.h b/Geo/gmshFace.h
index 59a368be3e6870111e8ef91913e45191f0088f1f..34440a14ab6dc7f2744ff2d557d341ff88404857 100644
--- a/Geo/gmshFace.h
+++ b/Geo/gmshFace.h
@@ -33,6 +33,7 @@ class gmshFace : public GFace {
void *getNativePtr() const { return s; }
virtual SPoint2 parFromPoint(const SPoint3 &) const;
virtual void resetMeshAttributes();
+ virtual bool checkTopology() {return true;};
};
#endif
diff --git a/Mesh/Generator.cpp b/Mesh/Generator.cpp
index 676cc757c7a4a8d346552b530799ef08986292f7..45db31cc540288042e48045a8b1ca11e774ac6ff 100644
--- a/Mesh/Generator.cpp
+++ b/Mesh/Generator.cpp
@@ -409,7 +409,12 @@ static void Mesh2D(GModel *m)
// and curve meshes) is global as it depends on a smooth normal
// field generated from the surface mesh of the source surfaces
if(!Mesh2DWithBoundaryLayers(m)){
- std::for_each(m->firstFace(), m->lastFace(), meshGFace());
+ //std::for_each(m->firstFace(), m->lastFace(), meshGFace());
+ std::set<GFace*> actualFaces;
+ for(GModel::fiter it = m->firstFace(); it != m->lastFace(); ++it){
+ actualFaces.insert(*it);
+ }
+ std::for_each(actualFaces.begin(), actualFaces.end(), meshGFace());
int nIter = 0;
while(1){
meshGFace mesher;
diff --git a/Mesh/meshGFace.cpp b/Mesh/meshGFace.cpp
index e877eff47a95f4d1568c25c913e323ea3c85ea9f..6d93676c53e2ec3d15b4ee926e4d726ab35b4a28 100644
--- a/Mesh/meshGFace.cpp
+++ b/Mesh/meshGFace.cpp
@@ -252,6 +252,8 @@ static bool meshGenerator(GFace *gf, int RECUR_ITER,
// replace edges by their compounds
if(gf->geomType() == GEntity::CompoundSurface){
+ bool correct = gf->checkTopology();
+ if (!correct) return true;
std::set<GEdge*> mySet;
std::list<GEdge*>::iterator it = edges.begin();
while(it != edges.end()){
diff --git a/Mesh/meshPartition.cpp b/Mesh/meshPartition.cpp
index 342e929221b70e765c79adc2a700775ce6bce147..afba861b7bd8e996e4221cb3782273ef3ba1553b 100644
--- a/Mesh/meshPartition.cpp
+++ b/Mesh/meshPartition.cpp
@@ -24,6 +24,9 @@
#include "partitionVertex.h"
#include "partitionEdge.h"
#include "partitionFace.h"
+#include "discreteEdge.h"
+#include "discreteFace.h"
+#include "GFaceCompound.h"
//--Prototype for Chaco interface
@@ -136,8 +139,8 @@ int PartitionMesh(GModel *const model, meshPartitionOptions &options)
model->recomputeMeshPartitions();
- // if (options.createPartitionBoundaries)
- CreatePartitionBoundaries (model);
+ if (options.createPartitionBoundaries)
+ CreatePartitionBoundaries (model);
Msg::Info("Partitioning complete");
Msg::StatusBar(1, false, "Mesh");
@@ -247,11 +250,18 @@ int PartitionGraph(Graph &graph, meshPartitionOptions &options)
metisOptions[2] = 1;
metisOptions[3] = options.refine_algorithm;
metisOptions[4] = 0;
- METIS_PartGraphKway
+ METIS_PartGraphKway
(&n, &graph.xadj[graph.section[iSec]],
&graph.adjncy[graph.section[iSec]], NULL, NULL, &wgtflag, &numflag,
&options.num_partitions, metisOptions, &edgeCut,
&graph.partition[graph.section[iSec]]);
+// printf("METIS with weights\n");
+// wgtflag = 2;
+// METIS_PartGraphKway
+// (&n, &graph.xadj[graph.section[iSec]],
+// &graph.adjncy[graph.section[iSec]], &graph.vwgts[graph.section[iSec]], NULL, &wgtflag, &numflag,
+// &options.num_partitions, metisOptions, &edgeCut,
+// &graph.partition[graph.section[iSec]]);
break;
}
}
@@ -724,7 +734,7 @@ void assignPartitionBoundary(GModel *model,
ppe = new partitionFace(model, -(int)pfaces.size()-1, v2);
pfaces.insert (ppe);
model->add(ppe);
- printf("created partitionFace %d (", ppe->tag());
+ printf("*** Created partitionFace %d (", ppe->tag());
for (unsigned int i = 0; i < v2.size(); ++i) printf("%d ", v2[i]);
printf(")\n");
}
@@ -765,7 +775,7 @@ void assignPartitionBoundary(GModel *model,
ppe = new partitionEdge(model, -(int)pedges.size()-1, 0, 0, v2);
pedges.insert (ppe);
model->add(ppe);
- printf("created partitionEdge %d (", ppe->tag());
+ printf("*** Create partitionEdge %d (", ppe->tag());
for (unsigned int i = 0; i < v2.size(); ++i) printf("%d ", v2[i]);
printf(")\n");
}
@@ -773,6 +783,72 @@ void assignPartitionBoundary(GModel *model,
ppe->lines.push_back(new MLine (me.getVertex(0),me.getVertex(1)));
}
+void splitBoundaryEdges(GModel *model, std::set<partitionEdge*, Less_partitionEdge> &newEdges)
+{
+
+ for (std::set<partitionEdge*, Less_partitionEdge>::iterator it = newEdges.begin() ; it != newEdges.end() ; ++it){
+
+ int nbSplit = 0;
+ partitionEdge *ge = *it;
+ std::list<MLine*> segments;
+ for (unsigned int i = 0; i < ge->lines.size(); i++){
+ segments.push_back(ge->lines[i]);
+ }
+
+ while (!segments.empty()) {
+ std::vector<MLine*> myLines;
+ std::list<MLine*>::iterator it = segments.begin();
+ MVertex *vB = (*it)->getVertex(0);
+ MVertex *vE = (*it)->getVertex(1);
+ myLines.push_back(*it);
+ segments.erase(it);
+ it++;
+ for (int i=0; i<2; i++) {
+ for (std::list<MLine*>::iterator it = segments.begin() ; it != segments.end(); ++it){
+ MVertex *v1 = (*it)->getVertex(0);
+ MVertex *v2 = (*it)->getVertex(1);
+ std::list<MLine*>::iterator itp;
+ if ( v1 == vE ){
+ myLines.push_back(*it);
+ itp = it;
+ it++;
+ segments.erase(itp);
+ vE = v2;
+ i = -1;
+ }
+ else if ( v2 == vE){
+ myLines.push_back(*it);
+ itp = it;
+ it++;
+ segments.erase(itp);
+ vE = v1;
+ i=-1;
+ }
+ if (it == segments.end()) break;
+ }
+ if (vB == vE) break;
+ if (segments.empty()) break;
+ MVertex *temp = vB;
+ vB = vE;
+ vE = temp;
+ }
+ if (nbSplit == 0 && segments.empty()) break;
+ int numEdge = model->maxEdgeNum() + 1;
+ discreteEdge *newGe = new discreteEdge(model, numEdge, 0, 0);
+ newGe->lines.insert(newGe->lines.end(), myLines.begin(), myLines.end());
+ model->add(newGe);
+ newGe->orderMLines(); //this creates also mesh_vertices
+
+ nbSplit++;
+ printf("*** split partitionEdge with tag =%d\n", numEdge);
+ }
+ if (nbSplit > 0) model->remove(ge);
+ }
+
+ return;
+
+}
+
void assignPartitionBoundary(GModel *model,
MVertex *ve,
std::set<partitionVertex*, Less_partitionVertex> &pvertices,
@@ -809,7 +885,7 @@ 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());
+ printf("*** created partitionVertex %d (", ppv->tag());
for (unsigned int i = 0; i < v2.size(); ++i) printf("%d ", v2[i]);
printf(")\n");
}
@@ -883,6 +959,7 @@ int CreatePartitionBoundaries(GModel *model)
}while (oper (e,it->first));
assignPartitionBoundary (model,e,pedges,voe,pfaces);
}
+ //splitBoundaryEdges(model,pedges);
}
}
@@ -922,6 +999,55 @@ int CreatePartitionBoundaries(GModel *model)
return 1;
}
+
+void CreateTopologyFromPartition(GModel *model, GFaceCompound *gf, int N)
+{
+
+ printf("---> CreateTopologyFromPartition for Compound Face %d \n", gf->tag());
+ // Compound is partitioned in N discrete faces
+ //--------------------------------------------
+ std::vector<discreteFace*> discreteFaces;
+ std::vector<std::set<MVertex*> > allNodes;
+ int numMax = model->maxFaceNum() + 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);
+ 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){
+ 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);
+ if (v0->onWhat()->dim() == 2) allNodes[part-1].insert(v0);
+ }
+ discreteFaces[part-1]->triangles.push_back(new MTriangle(e->getVertex(0),e->getVertex(1),e->getVertex(2)));
+ }
+ }
+
+ for( int i =0; i < N; i++){
+ GFace *face = model->getFaceByTag(numMax+i);
+ for (std::set<MVertex*>::iterator it = allNodes[i].begin(); it != allNodes[i].end(); it++){
+ face->mesh_vertices.push_back(*it);
+ }
+ }
+
+ //remove the discrete face that is not topologically correct
+ //for(it = _compound.begin() ; it != _compound.end() ; ++it){
+ // model->remove(*it);
+ //}
+
+ return;
+
+}
/*******************************************************************************
*
diff --git a/Mesh/meshPartition.h b/Mesh/meshPartition.h
index b9d40cdea91736285cfc11fa5a7c12df7ebb70de..2b8728f6b4cc6320a194c4b452d99fce2fe6c932 100644
--- a/Mesh/meshPartition.h
+++ b/Mesh/meshPartition.h
@@ -7,6 +7,8 @@
#define _PARTITION_H_
#include <vector>
+#include "partitionEdge.h"
+#include "GFaceCompound.h"
struct meshPartitionOptions;
struct BoElemGr;
@@ -26,5 +28,7 @@ int MakeGraph(GModel *const model, Graph &graph,
int PartitionGraph(Graph &graph, meshPartitionOptions &options);
int PartitionMesh(GModel *const model, meshPartitionOptions &options);
int CreatePartitionBoundaries (GModel *model);
+void splitBoundaryEdges(GModel *model, std::set<partitionEdge*, Less_partitionEdge> &newEdges);
+void CreateTopologyFromPartition(GModel *model, GFaceCompound *gf, int num_parts);
#endif
diff --git a/Mesh/meshPartitionObjects.h b/Mesh/meshPartitionObjects.h
index 96680252d59f696bac6ff79acd462e59da9cc127..309ed488c59c12816da4a1ef57c018af16136541 100644
--- a/Mesh/meshPartitionObjects.h
+++ b/Mesh/meshPartitionObjects.h
@@ -29,6 +29,7 @@ class GrVertex
public:
const int index; // This is the creation index, *not* the
// index in 'adjncy'
+ int dualWeight;
private:
unsigned short size;
unsigned short sizeC; // Complete size (all possible 'grEdge')
@@ -86,6 +87,8 @@ class Graph
std::vector<int> adjncy; // Connectivity between graph vertex
// xadj[i] and its neighbour graph
// vertices.
+ std::vector<int> vwgts; // Weights assigned for each
+ // vertex
std::vector<int> section; // For separate partitioning of
// different parts of the mesh
std::vector<int> partition; // The partitions output from the
@@ -128,6 +131,7 @@ class Graph
totalGrVert = _totalGrVert;
xadj.resize(_totalGrVert + 1);
adjncy.reserve(2*totalGrEdge);
+ vwgts.resize(_totalGrVert);
partition.resize(_totalGrVert);
element.resize(_totalGrVert);
c2w = new int[_totalGrVert];
@@ -137,6 +141,7 @@ class Graph
{
const int i = numGrVert++;
xadj[i] = adjncy.size();
+ vwgts[i-1]=(int)(1.0);
grVertMapIt->second.write(adjncy);
element[i] = grVertMapIt->first;
// Translated vertex numbers start from 1
@@ -154,6 +159,7 @@ class Graph
Msg::Warning("Internal error - Graph vertices are missing");
}
xadj[numGrVert] = adjncy.size();
+ vwgts[numGrVert-1]=(int)(1.0);
const int nAdj = adjncy.size();
for(int i = 0; i != nAdj; ++i) adjncy[i] = c2w[adjncy[i]];
delete[] c2w;
diff --git a/Solver/convexCombinationTerm.h b/Solver/convexCombinationTerm.h
index 1f77c8137ac72423a8ad520c9cafc7c2e8fca7dd..dcaef8fb3b8f7546895c5b60181a8a1bd6361c39 100644
--- a/Solver/convexCombinationTerm.h
+++ b/Solver/convexCombinationTerm.h
@@ -6,20 +6,22 @@
#ifndef _CONVEX_COMBINATION_TERM_H_
#define _CONVEX_COMBINATION_TERM_H_
+#include <assert.h>
#include "femTerm.h"
#include "simpleFunction.h"
#include "Gmsh.h"
#include "GModel.h"
#include "SElement.h"
#include "fullMatrix.h"
+#include "Numeric.h"
-class convexCombinationTerm : public femTerm<double, double> {
+class convexCombinationTerm : public femTerm<double,double> {
protected:
- const simpleFunction<double> *_diffusivity;
+ const simpleFunction<double> *_k;
const int _iField;
public:
- convexCombinationTerm(GModel *gm, int iField, simpleFunction<double> *diffusivity)
- : femTerm<double, double>(gm), _iField(iField), _diffusivity(diffusivity) {}
+ convexCombinationTerm(GModel *gm, int iField, simpleFunction<double> *k)
+ : femTerm<double,double>(gm), _iField(iField), _k(k) {}
virtual int sizeOfR(SElement *se) const
{
return se->getMeshElement()->getNumVertices();
@@ -30,16 +32,16 @@ class convexCombinationTerm : public femTerm<double, double> {
}
Dof getLocalDofR(SElement *se, int iRow) const
{
- return Dof(se->getMeshElement()->getVertex(iRow)->getNum(), _iField);
+ return Dof(se->getMeshElement()->getVertex(iRow)->getNum(),
+ Dof::createTypeWithTwoInts(0, _iField));
}
- virtual void elementMatrix(SElement *e, fullMatrix<double> &m) const
+ virtual void elementMatrix(SElement *se, fullMatrix<double> &m) const
{
+
+ MElement *e = se->getMeshElement();
m.set_all(0.);
- int nbNodes = e->getMeshElement()->getNumVertices();
- //double x = 0.3*(e->getVertex(0)->x()+e->getVertex(1)->x()+e->getVertex(2)->x());
- //double y = 0.3*(e->getVertex(0)->y()+e->getVertex(1)->y()+e->getVertex(2)->y());
- //double z = 0.3*(e->getVertex(0)->z()+e->getVertex(1)->z()+e->getVertex(2)->z());
- const double _diff = 1.0;
+ const int nbNodes = e->getNumVertices();
+ const double _diff = 1.0;
for (int j = 0; j < nbNodes; j++){
for (int k = 0; k < nbNodes; k++){
m(j,k) = -1.*_diff;
@@ -47,6 +49,8 @@ class convexCombinationTerm : public femTerm<double, double> {
m(j,j) = (nbNodes - 1) * _diff;
}
}
+
+
};
#endif
diff --git a/Solver/helmholtzTerm.h b/Solver/helmholtzTerm.h
index 7eb14217710f5ea672aaf572fb10f437e965a852..c120bbb2adff8e596346f80deaa18d9e13715c44 100644
--- a/Solver/helmholtzTerm.h
+++ b/Solver/helmholtzTerm.h
@@ -15,7 +15,7 @@
#include "fullMatrix.h"
#include "Numeric.h"
-// \nabla \cdot k \nabla U + a U
+// \nabla \cdot k \nabla U - a U
template<class scalar>
class helmholtzTerm : public femTerm<scalar, scalar> {
protected:
diff --git a/benchmarks/2d/square.geo b/benchmarks/2d/square.geo
index 34634b39110208370f6c4029e643a3cf52fc3679..6475d2fad2650de9072f870f760f2839d37d3457 100644
--- a/benchmarks/2d/square.geo
+++ b/benchmarks/2d/square.geo
@@ -11,5 +11,5 @@ Line(4) = {1, 2};
Line Loop(5) = {1, 2, 3, 4};
Plane Surface(10) = {5};
-Compound Surface(11)={10} Boundary{{}};
+Compound Surface(11)={10};
diff --git a/benchmarks/3d/Torus2.geo b/benchmarks/3d/Torus2.geo
index 0460cb7ea84b609456c263cd3fff2e5cae7a480b..245f082097f7e64738fb512962545044168a83d3 100644
--- a/benchmarks/3d/Torus2.geo
+++ b/benchmarks/3d/Torus2.geo
@@ -14,4 +14,7 @@ Extrude Surface {6, {0,1,0}, {-2,0,0}, 2*Pi/3};
Extrude Surface {28, {0,1,0}, {-2,0,0}, 2*Pi/3};
Extrude Surface {50, {0,1,0}, {-2,0,0}, 2*Pi/3};
Surface Loop(72) = {45,23,67,71,49,27,15,59,37,41,19,63};
-Volume(73) = {72};
+
+//Compound Surface(100)={45,23,67,71,49,27,15,59,37,41,19,63};
+
+//Volume(73) = {72};
diff --git a/benchmarks/3d/Torus_GEO.geo b/benchmarks/3d/Torus_GEO.geo
index c0dc7f1354e026c6701d2a155880563e5feb4fd9..bdda121a002a85fac8b4d8aec41bdd98b5971184 100644
--- a/benchmarks/3d/Torus_GEO.geo
+++ b/benchmarks/3d/Torus_GEO.geo
@@ -1,4 +1,4 @@
-Mesh.CharacteristicLengthFactor=1.0;
+Mesh.CharacteristicLengthFactor=0.4;
Merge "Torus2_CLASS.msh";
CreateTopology;
@@ -6,5 +6,6 @@ CreateTopology;
Compound Line(100)={1};
Compound Line(200)={2};
-Compound Surface(1000)={150} Boundary {{}};
-Compound Surface(2000)={250} Boundary {{}};
+Compound Surface(1000)={150};
+Compound Surface(2000)={250};
+
diff --git a/benchmarks/stl/PelvisARTHUR_CLASS_GEO.geo b/benchmarks/stl/PelvisARTHUR_CLASS_GEO.geo
index 3e288e9ac468e3b813bbffc8e4b4cea5e8a73912..437b2d692850cb7cbed47a5cbdf4ec2dfea48659 100644
--- a/benchmarks/stl/PelvisARTHUR_CLASS_GEO.geo
+++ b/benchmarks/stl/PelvisARTHUR_CLASS_GEO.geo
@@ -1,13 +1,13 @@
-Mesh.CharacteristicLengthFactor=0.1;
+Mesh.CharacteristicLengthFactor=0.2;
-Merge "PelvisARTHUR_CLASS.msh";
+Merge "PelvisKO.msh"; //ARTHUR_CLASS.msh";
CreateTopology;
Surface Loop(1002)={2,3};
Volume(1003)={1002};
Compound Line(10)={4};
-Compound Line(20)={5};
+//Compound Line(20)={5};
Compound Surface(200)={2};
Compound Surface(400)={3};
diff --git a/benchmarks/stl/implant_CLASS_GEO.geo b/benchmarks/stl/implant_CLASS_GEO.geo
index 8636edf12c00bf61ac247e2a83c1896f58779c36..fcfa5a5a1766c3fb188f9f4cfe1dc5a8331f8229 100644
--- a/benchmarks/stl/implant_CLASS_GEO.geo
+++ b/benchmarks/stl/implant_CLASS_GEO.geo
@@ -1,4 +1,4 @@
-Mesh.CharacteristicLengthFactor=0.5;
+Mesh.CharacteristicLengthFactor=0.1;
Merge "implant_CLASS.msh";
CreateTopology;
@@ -7,7 +7,6 @@ Compound Line(20)={3};
Compound Surface(100)={2} Boundary {{}};
Compound Surface(200)={3} Boundary {{}};
-
//Compound Line(20)={6};
//Compound Line(30)={7};
//Compound Line(40)={8};