From c65d59fbcfa5c9cf4657b20aeb9d7a07b00b21bb Mon Sep 17 00:00:00 2001
From: Emilie Marchandise <emilie.marchandise@uclouvain.be>
Date: Wed, 18 Nov 2009 19:40:43 +0000
Subject: [PATCH]
---
Geo/GFaceCompound.cpp | 176 ++++++++++++++++++++++-----------
Geo/GFaceCompound.h | 9 +-
Geo/GModel.cpp | 2 +-
Geo/SOrientedBoundingBox.cpp | 5 +
Geo/SOrientedBoundingBox.h | 1 +
Mesh/meshPartition.cpp | 17 +---
Mesh/multiscalePartition.cpp | 135 ++++++++++++++++++++-----
Solver/convexCombinationTerm.h | 1 +
8 files changed, 249 insertions(+), 97 deletions(-)
diff --git a/Geo/GFaceCompound.cpp b/Geo/GFaceCompound.cpp
index 6afed65789..aae88152ff 100644
--- a/Geo/GFaceCompound.cpp
+++ b/Geo/GFaceCompound.cpp
@@ -188,6 +188,59 @@ static void myPolygon(std::vector<MElement*> &vTri, std::vector<MVertex*> &vPoly
}
+void GFaceCompound::fillNeumannBCS() const
+{
+
+ //close neuman bcs
+ for(std::list<std::list<GEdge*> >::const_iterator iloop = _interior_loops.begin();
+ iloop != _interior_loops.end(); iloop++){
+ std::list<GEdge*> loop = *iloop;
+ if (loop != _U0 ){
+ //--- center of Neumann interior loop
+ int nb = 0;
+ double x=0.;
+ double y=0.;
+ double z=0.;
+ for (std::list<GEdge*>::iterator ite = loop.begin(); ite != loop.end(); ite++){
+ for (int k= 0; k< (*ite)->getNumMeshElements(); k++){
+ MVertex *v0 = (*ite)->getMeshElement(k)->getVertex(0);
+ MVertex *v1 = (*ite)->getMeshElement(k)->getVertex(1);
+ x += .5*(v0->x() + v1->x());
+ y += .5*(v0->y() + v1->y());
+ z += .5*(v0->z() + v1->z());
+ nb++;
+ }
+ }
+ x/=nb; y/=nb; z/=nb;
+ MVertex *c = new MVertex(x, y, z);
+
+ //--- create new triangles
+ for (std::list<GEdge*>::iterator ite = loop.begin(); ite != loop.end(); ite++){
+ for (int i= 0; i< (*ite)->getNumMeshElements(); i++){
+ MVertex *v0 = (*ite)->getMeshElement(i)->getVertex(0);
+ MVertex *v1 = (*ite)->getMeshElement(i)->getVertex(1);
+ MTriangle *myTri = new MTriangle(v0,v1, c);
+ fillTris.push_back(myTri);
+ }
+ }
+ }
+ }
+
+// FILE * ftri = fopen("fillTris.pos","w");
+// fprintf(ftri,"View \"\"{\n");
+// for (std::list<MTriangle*>::iterator it2 = fillTris.begin(); it2 !=fillTris.end(); it2++ ){
+// MTriangle *t = (*it2);
+// fprintf(ftri,"ST(%g,%g,%g,%g,%g,%g,%g,%g,%g){%g,%g,%g};\n",
+// t->getVertex(0)->x(), t->getVertex(0)->y(), t->getVertex(0)->z(),
+// t->getVertex(1)->x(), t->getVertex(1)->y(), t->getVertex(1)->z(),
+// t->getVertex(2)->x(), t->getVertex(2)->y(), t->getVertex(2)->z(),
+// 1., 1., 1.);
+// }
+// fprintf(ftri,"};\n");
+// fclose(ftri);
+
+}
+
bool GFaceCompound::trivial() const
{
if(_compound.size() == 1 &&
@@ -255,7 +308,7 @@ bool GFaceCompound::checkOrientation(int iter) const
//Only check orientation for stl files (1 patch)
// if(_compound.size() > 1.0) return true;
- return true;
+ //return true;
std::list<GFace*>::const_iterator it = _compound.begin();
double a_old = 0, a_new;
@@ -370,34 +423,6 @@ void GFaceCompound::one2OneMap() const
SPoint3 p_cg(u_cg,v_cg,0);
coordinates[v] = p_cg;
-// 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");
-
-// for(int i=0; i< vTri.size(); i++){
-// MTriangle *t = (MTriangle*) vTri[i];
-// SPoint3 v1 = coordinates[t->getVertex(0)];
-// SPoint3 v2 = coordinates[t->getVertex(1)];
-// SPoint3 v3 = coordinates[t->getVertex(2)];
-
-// printf("//////////////////// \n " );
-// double p1[2] = {v1[0],v1[1]};
-// double p2[2] = {v2[0],v2[1]};
-// double p3[2] = {v3[0],v3[1]};
-// double a_new = robustPredicates::orient2d(p1, p2, p3);
-// MVertex *e1=t->getVertex(0); MVertex *e2=t->getVertex(1); MVertex *e3=t->getVertex(2);
-// printf("Point(%d)={%g, %g, 0}; \n ",e1->getNum(), v1[0],v1[1] );
-// printf("Point(%d)={%g, %g, 0}; \n ",e2->getNum(), v2[0],v2[1] );
-// printf("Point(%d)={%g, %g, 0}; \n ",e3->getNum(), v3[0],v3[1] );
-// printf("Line(%d)={%d,%d}; \n", e1->getNum()+e2->getNum() , e1->getNum(), e2->getNum());
-// printf("Line(%d)={%d,%d}; \n", e2->getNum()+e3->getNum() , e2->getNum(), e3->getNum());
-// printf("Line(%d)={%d,%d}; \n", e3->getNum()+e1->getNum() , e3->getNum(), e1->getNum());
-// printf("Surface(%d)={%d,%d, %d}; \n", e3->getNum()+e1->getNum() + e2->getNum(), e3->getNum(), e1->getNum()+e2->getNum() ,
-// e2->getNum()+e3->getNum() , e3->getNum()+e1->getNum() );
-// printf("//Area=%g \n", a_new);
-// }
-
}
}
#endif
@@ -418,7 +443,7 @@ bool GFaceCompound::parametrize() const
//-----------------
if (_mapping == HARMONIC){
Msg::Info("Parametrizing surface %d with 'harmonic map'", tag());
- parametrize(ITERU,HARMONIC);
+ parametrize(ITERU,HARMONIC);
parametrize(ITERV,HARMONIC);
}
//Conformal map parametrization
@@ -441,8 +466,11 @@ bool GFaceCompound::parametrize() const
if (!checkOrientation(0)){
Msg::Info("Parametrization failed using standard techniques : moving to convex combination");
coordinates.clear();
+ Octree_Delete(oct);
+ fillNeumannBCS();
parametrize(ITERU,CONVEXCOMBINATION);
parametrize(ITERV,CONVEXCOMBINATION);
+ buildOct();
}
computeNormals();
@@ -464,7 +492,7 @@ void GFaceCompound::getBoundingEdges()
l_edges.clear();
if(_U0.size()){
- //in case the bounding edges are explicitely given
+ //--- 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);
@@ -480,40 +508,64 @@ void GFaceCompound::getBoundingEdges()
while(!_unique.empty()) computeALoop(_unique,loop);
}
else{
- //in case the bounding edges are NOT explicitely given
+ //--- in case the bounding edges are NOT explicitely given
std::set<GEdge*>::iterator itf = _unique.begin();
for( ; itf != _unique.end(); ++itf){
- //printf("for Compound face %d add U0 bounding edge %d \n", tag(), (*itf)->tag());
l_edges.push_back(*itf);
(*itf)->addFace(this);
}
- computeALoop(_unique, _U0);
- while(!_unique.empty()) computeALoop(_unique,_U1);
+ std::list<GEdge*> loop;
+ computeALoop(_unique,loop); //_U0);
+ while(!_unique.empty()) computeALoop(_unique, loop); //_U1);
+
+ //assign Derichlet BC (_U0) to bound with largest size
+ double maxSize = 0.0;
+ for(std::list<std::list<GEdge*> >::iterator it = _interior_loops.begin();
+ it != _interior_loops.end(); it++){
+ double size = getSizeBB(*it);
+ if (size > maxSize) {
+ _U0 = *it;
+ maxSize = size;
+ }
+ }
+
}
return;
}
-SBoundingBox3d GFaceCompound::bound_U0() const
+
+double GFaceCompound::getSizeBB(const std::list<GEdge* > &elist) const
+{
+
+ SOrientedBoundingBox obboxD = obb_boundEdges(elist);
+ double D = obboxD.getMaxSize();
+
+ //SOrientedBoundingBox bboxD = boundEdges(elist);
+ //double D = norm(SVector3(bboxD.max(), bboxD.min()));
+
+ return D;
+
+}
+SBoundingBox3d GFaceCompound::boundEdges(const std::list<GEdge* > &elist) const
{
SBoundingBox3d res;
- std::list<GEdge*>::const_iterator it = _U0.begin();
- for(; it != _U0.end(); it++)
+ std::list<GEdge*>::const_iterator it = elist.begin();
+ for(; it != elist.end(); it++)
res += (*it)->bounds();
return res;
}
-SOrientedBoundingBox GFaceCompound::obb_bound_U0() const
+SOrientedBoundingBox GFaceCompound::obb_boundEdges(const std::list<GEdge* > &elist) const
{
SOrientedBoundingBox res;
std::vector<SPoint3> vertices;
-
-std::list<GEdge*>::const_iterator it = _U0.begin();
-for(; it != _U0.end(); it++) {
-
+ std::list<GEdge*>::const_iterator it = elist.begin();
+ for(; it != elist.end(); it++) {
+
if((*it)->getNumMeshVertices() > 0) {
int N = (*it)->getNumMeshVertices();
for (int i = 0; i < N; i++) {
@@ -921,7 +973,16 @@ void GFaceCompound::parametrize(iterationStep step, typeOfMapping tom) const
myAssembler.numberVertex(t->getVertex(1), 0, 1);
myAssembler.numberVertex(t->getVertex(2), 0, 1);
}
- }
+ }
+ if (tom == CONVEXCOMBINATION){
+ for (std::list<MTriangle*>::iterator it2 = fillTris.begin(); it2 !=fillTris.end(); it2++ ){
+ MTriangle *t = (*it2);
+ myAssembler.numberVertex(t->getVertex(0), 0, 1);
+ myAssembler.numberVertex(t->getVertex(1), 0, 1);
+ myAssembler.numberVertex(t->getVertex(2), 0, 1);
+ }
+ }
+
Msg::Debug("Creating term %d dofs numbered %d fixed",
myAssembler.sizeOfR(), myAssembler.sizeOfF());
@@ -937,6 +998,10 @@ void GFaceCompound::parametrize(iterationStep step, typeOfMapping tom) const
laplace.addToMatrix(myAssembler, &se);
}
}
+ for (std::list<MTriangle*>::iterator it2 = fillTris.begin(); it2 !=fillTris.end(); it2++ ){
+ SElement se((*it2));
+ laplace.addToMatrix(myAssembler, &se);
+ }
}
else {
laplaceTerm laplace(model(), 1, &ONE);
@@ -1046,8 +1111,6 @@ void GFaceCompound::parametrize_conformal() const
double value2 = myAssembler.getDofValue(v,0,2);
coordinates[v] = SPoint3(value1,value2,0.0);
}
- // printStuff();
- //exit(1);
_lsys->clear();
@@ -1360,7 +1423,7 @@ void GFaceCompound::getTriangle(double u, double v,
void GFaceCompound::buildOct() const
{
printStuff();
-
+
SBoundingBox3d bb;
int count = 0;
std::list<GFace*>::const_iterator it = _compound.begin();
@@ -1424,8 +1487,8 @@ void GFaceCompound::buildOct() const
bool GFaceCompound::checkTopology() const
{
// FIXME!!! I think those things are wrong with cross-patch reparametrization
- //if ((*(_compound.begin()))->geomType() != GEntity::DiscreteSurface)return true;
-
+ //if ((*(_compound.begin()))->geomType() != GEntity::DiscreteSurface)return true;
+
bool correctTopo = true;
int Nb = _interior_loops.size();
@@ -1434,7 +1497,7 @@ bool GFaceCompound::checkTopology() const
if (!correctTopo){
Msg::Warning("Wrong topology: Genus=%d and N boundary=%d", G, Nb);
- nbSplit = std::max(G+1, 2);
+ nbSplit = G+2; //std::max(G+2, 2);
Msg::Info("Split surface %d in %d parts with Mesh partitioner", tag(), nbSplit);
}
else
@@ -1452,7 +1515,7 @@ bool GFaceCompound::checkAspectRatio() const
bool paramOK = true;
if(allNodes.empty()) buildAllNodes();
- double limit = 1.e17 ;
+ double limit = 1.e15 ;
double areaMax = 0.0;
int nb = 0;
std::list<GFace*>::const_iterator it = _compound.begin();
@@ -1483,12 +1546,11 @@ bool GFaceCompound::checkAspectRatio() const
if (areaMax > limit && nb > 10) {
Msg::Warning("Geometrical aspect ratio too high (1/area_2D=%g)", areaMax);
SBoundingBox3d bboxH = bounds();
- SOrientedBoundingBox obboxD = obb_bound_U0();
- double H = norm(SVector3(bboxH.max(), bboxH.min()));
- double D = obboxD.getMaxSize();
- int split = (int)ceil(.3*H/D);
+ double H = norm(SVector3(bboxH.max(), bboxH.min()));
+ double D = getSizeBB(_U0);
+ int split = (int)ceil(.7*H/D);
nbSplit = std::max(split,2);
- printf("H=%g, D=%g split=%d nbSplit=%d \n", H, D, split, nbSplit);
+ //printf("H=%g, D=%g split=%d nbSplit=%d \n", H, D, split, nbSplit);
Msg::Info("Partition geometry in N=%d parts", nbSplit);
paramOK = false;
}
@@ -1594,7 +1656,7 @@ int GFaceCompound::genusGeom() const
void GFaceCompound::printStuff() const
{
std::list<GFace*>::const_iterator it = _compound.begin();
-
+
char name0[256], name1[256], name2[256], name3[256];
char name4[256], name5[256], name6[256];
sprintf(name0, "UVAREA-%d.pos", (*it)->tag());
diff --git a/Geo/GFaceCompound.h b/Geo/GFaceCompound.h
index 73c44debe8..7efa430cb5 100644
--- a/Geo/GFaceCompound.h
+++ b/Geo/GFaceCompound.h
@@ -66,6 +66,7 @@ class GFaceCompound : public GFace {
mutable bool mapv2Tri;
mutable std::map<MVertex*, SPoint3> coordinates;
mutable std::map<MVertex*, SVector3> _normals;
+ mutable std::list<MTriangle*> fillTris;
void buildOct() const ;
void buildAllNodes() const;
void parametrize(iterationStep, typeOfMapping) const;
@@ -85,9 +86,11 @@ class GFaceCompound : public GFace {
void printStuff() const;
bool trivial() const ;
linearSystem <double> *_lsys;
- SBoundingBox3d bound_U0() const;
- SOrientedBoundingBox obb_bound_U0() const;
- public:
+ double getSizeBB(const std::list<GEdge* > &elist) const;
+ SBoundingBox3d boundEdges(const std::list<GEdge* > &elist) const;
+ SOrientedBoundingBox obb_boundEdges(const std::list<GEdge* > &elist) const;
+ void fillNeumannBCS() const;
+ public:
GFaceCompound(GModel *m, int tag, std::list<GFace*> &compound,
std::list<GEdge*> &U0, std::list<GEdge*> &U1,
std::list<GEdge*> &V0, std::list<GEdge*> &V1,
diff --git a/Geo/GModel.cpp b/Geo/GModel.cpp
index eee25937eb..aece6243d6 100644
--- a/Geo/GModel.cpp
+++ b/Geo/GModel.cpp
@@ -1230,7 +1230,7 @@ void GModel::createTopologyFromFaces(std::vector<discreteFace*> &discFaces)
int numE = maxEdgeNum()+1;
discreteEdge *e = new discreteEdge(this, numE, 0, 0);
- printf("*** Created discreteEdge %d \n", numE);
+ //printf("*** Created discreteEdge %d \n", numE);
add(e);
discEdges.push_back(e);
diff --git a/Geo/SOrientedBoundingBox.cpp b/Geo/SOrientedBoundingBox.cpp
index 7a9dad595e..128ae0637f 100644
--- a/Geo/SOrientedBoundingBox.cpp
+++ b/Geo/SOrientedBoundingBox.cpp
@@ -152,6 +152,11 @@ double SOrientedBoundingBox::getMaxSize()
return (std::max(size[0], std::max(size[1], size[2])));
}
+double SOrientedBoundingBox::getMinSize()
+{
+ return (std::min(size[0], std::min(size[1], size[2])));
+}
+
SVector3 SOrientedBoundingBox::getAxis(int axis)
{
SVector3 ret;
diff --git a/Geo/SOrientedBoundingBox.h b/Geo/SOrientedBoundingBox.h
index e2fd50b979..2d9265772e 100644
--- a/Geo/SOrientedBoundingBox.h
+++ b/Geo/SOrientedBoundingBox.h
@@ -64,6 +64,7 @@ class SOrientedBoundingBox {
double getCenterZ(){ return center[2]; }
SVector3 getSize(){ return size; }
double getMaxSize();
+ double getMinSize();
// valid values for axis are 0 (X-axis), 1 (Y-axis) or 2 (Z-axis)
SVector3 getAxis(int axis);
diff --git a/Mesh/meshPartition.cpp b/Mesh/meshPartition.cpp
index ddcd14e326..b9707b7b8f 100644
--- a/Mesh/meshPartition.cpp
+++ b/Mesh/meshPartition.cpp
@@ -114,6 +114,8 @@ bool PartitionZeroGenus(std::list<GFace*> &cFaces, int &nbParts){
options.mesh_dims[0] = nbParts;
}
+// PartitionMeshFace(cFaces, options);
+
std::vector<MElement *> elements;
for (std::list<GFace*>::iterator it = cFaces.begin(); it != cFaces.end(); it++)
for(unsigned int j = 0; j < (*it)->getNumMeshElements(); j++)
@@ -124,19 +126,6 @@ bool PartitionZeroGenus(std::list<GFace*> &cFaces, int &nbParts){
}
-// bool PartitionZeroGenus(std::list<GFace*> &cFaces, int &nbParts){
-// meshPartitionOptions options;
-// options = CTX::instance()->partitionOptions;
-// options.num_partitions = nbParts;
-// options.partitioner = 1;//1 CHACO //2 METIS
-// if ( options.partitioner == 1){
-// options.global_method = 2;// 1 Multilevel-KL 2 Spectral
-// options.mesh_dims[0] = nbParts;
-// }
-// PartitionMeshFace(cFaces, options);
-// return true;
-// }
-
int PartitionMeshElements( std::vector<MElement*> &elements, meshPartitionOptions &options){
GModel *tmp_model = new GModel();
@@ -1097,7 +1086,7 @@ void createPartitionFaces(GModel *model, GFaceCompound *gf, int N,
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);
+ //printf("*** Created discreteFace %d \n", numMax+i);
discreteFace *face = new discreteFace(model, numMax+i);
discreteFaces.push_back(face);
model->add(face);//delete this
diff --git a/Mesh/multiscalePartition.cpp b/Mesh/multiscalePartition.cpp
index a772b78acc..84d282eaa3 100644
--- a/Mesh/multiscalePartition.cpp
+++ b/Mesh/multiscalePartition.cpp
@@ -2,12 +2,53 @@
#include "GmshConfig.h"
#include "GmshDefines.h"
#include "meshPartition.h"
+#include "MEdge.h"
#include "MElement.h"
-static bool zeroGenus (std::vector<MElement *> &elements){
+static void recur_connect (MVertex *v,
+ std::multimap<MVertex*,MEdge> &v2e,
+ std::set<MEdge,Less_Edge> &group,
+ std::set<MVertex*> &touched){
- //TODO
- return true;
+ if (touched.find(v) != touched.end())return;
+
+ touched.insert(v);
+ for (std::multimap <MVertex*,MEdge>::iterator it = v2e.lower_bound(v);
+ it != v2e.upper_bound(v) ; ++it){
+ group.insert(it->second);
+ for (int i=0;i<it->second.getNumVertices();++i){
+ recur_connect (it->second.getVertex(i),v2e,group,touched);
+ }
+ }
+
+}
+
+static int connected_bounds (std::vector<MEdge> &edges)
+{
+
+ std::vector<std::vector<MEdge> > regions;
+
+ std::multimap<MVertex*,MEdge> v2e;
+ for (int i=0;i<edges.size();++i){
+ for (int j=0;j<edges[i].getNumVertices();j++){
+ v2e.insert(std::make_pair(edges[i].getVertex(j),edges[i]));
+ }
+ }
+ while (!v2e.empty()){
+ std::set<MEdge, Less_Edge> group;
+ std::set<MVertex*> touched;
+ recur_connect (v2e.begin()->first,v2e,group,touched);
+ std::vector<MEdge> temp;
+ temp.insert(temp.begin(), group.begin(), group.end());
+ regions.push_back(temp);
+ for ( std::set<MVertex*>::iterator it = touched.begin() ; it != touched.end();++it)
+ v2e.erase(*it);
+ }
+
+ return regions.size();
+}
+
+static int getGenus (std::vector<MElement *> &elements){
//We suppose MElements are simply connected
@@ -28,19 +69,51 @@ static bool zeroGenus (std::vector<MElement *> &elements){
int poincare = vs.size() - es.size() + N;
//compute connected boundaries
- //....
int nbBounds = 0;
- //....
-
+ std::vector<MEdge> bEdges;
+ for(unsigned int i = 0; i < elements.size(); i++){
+ for(unsigned int j = 0; j < elements[i]->getNumEdges(); j++){
+ MEdge me = elements[i]->getEdge(j);
+ if(std::find(bEdges.begin(), bEdges.end(), me) == bEdges.end())
+ bEdges.push_back(me);
+ else
+ bEdges.erase(std::find(bEdges.begin(), bEdges.end(),me));
+ }
+ }
+ nbBounds = connected_bounds(bEdges);
int genus = (int)(-poincare + 2 - nbBounds)/2;
+
+ //printf("************** partition has %d boundaries and genus =%d \n", nbBounds, genus);
- if (genus == 0)
- return true;
- else
- return false;
+ return genus;
}
+// static int getGeomAspectRatio (std::vector<MElement *> &elements){
+
+// std::set<MVertex*> vs;
+// for(unsigned int i = 0; i < elements.size(); i++){
+// MElement *e = elements[i];
+// for(int j = 0; j < e->getNumVertices(); j++){
+// vs.insert(e->getVertex(j));
+// }
+// }
+
+// std::vector<SPoint3> vertices;
+// for (std::set<MVertex* >::iterator it = vs.begin(); it != vs.end(); it++){
+// SPoint3 pt((*it)->x(),(*it)->y(), (*it)->z());
+// vertices.push_back(pt);
+// }
+
+// SOrientedBoundingBox obbox = SOrientedBoundingBox::buildOBB(vertices);
+// double eta = obbox.getMaxSize()/obbox.getMinSize();
+
+// printf("aspect ratio = %g (max=%g, min=%g)\n", eta, obbox.getMaxSize(), obbox.getMinSize() );
+
+// return eta;
+
+// }
+
static void partitionRegions (std::vector<MElement*> &elements,
std::vector<std::vector<MElement*> > ®ions){
@@ -89,16 +162,21 @@ void multiscalePartition::partition(partitionLevel & level){
levels.push_back(nextLevel);
- if (!zeroGenus(regions[i])){
- Msg::Info("Multiscale partition, level %d region %d not ZERO-GENUS",
- nextLevel->recur,nextLevel->region);
- partition(*nextLevel);
- }
- else {
- Msg::Info("Multiscale Partition, level %d, region %d is ZERO-GENUS",
- nextLevel->recur,nextLevel->region);
- }
-
+ int genus = getGenus(regions[i]);
+ if (genus < 0) {
+ Msg::Error("Genus partition is negative G=%d!", genus);
+ exit(1);
+ }
+ if (genus != 0 ){
+ Msg::Info("Multiscale partition, level %d region %d is %d-GENUS -> partition",
+ nextLevel->recur,nextLevel->region, genus);
+ partition(*nextLevel);
+ }
+ else {
+ Msg::Info("Multiscale Partition, level %d, region %d is ZERO-GENUS",
+ nextLevel->recur,nextLevel->region);
+ }
+
}
#endif
@@ -106,8 +184,21 @@ void multiscalePartition::partition(partitionLevel & level){
int multiscalePartition::assembleAllPartitions(){
- int nbParts = options.num_partitions;;
+ int nbParts = 1;
+ //int nbLevel = options.num_partitions;
+
+ for (int i = 0; i< levels.size(); i++){
+ partitionLevel *iLevel = levels[i];
+ if(iLevel->elements.size() > 0){
+ for (int j = 0; j < iLevel->elements.size(); j++){
+ MElement *e = iLevel->elements[j];
+ int part = e->getPartition();
+ e->setPartition(nbParts);
+ }
+ nbParts++;
+ }
+ }
- return nbParts;
+ return nbParts-1;
}
diff --git a/Solver/convexCombinationTerm.h b/Solver/convexCombinationTerm.h
index 2e3e068f5b..d27de599cc 100644
--- a/Solver/convexCombinationTerm.h
+++ b/Solver/convexCombinationTerm.h
@@ -1,3 +1,4 @@
+
// Gmsh - Copyright (C) 1997-2009 C. Geuzaine, J.-F. Remacle
//
// See the LICENSE.txt file for license information. Please report all
--
GitLab