diff --git a/Geo/GFaceCompound.cpp b/Geo/GFaceCompound.cpp
index aae88152ff08426b16e16a3641588e454414fe50..305f283d3e462f4d6fd21e6a44a2bf3f8dbd9cb7 100644
--- a/Geo/GFaceCompound.cpp
+++ b/Geo/GFaceCompound.cpp
@@ -467,7 +467,7 @@ bool GFaceCompound::parametrize() const
Msg::Info("Parametrization failed using standard techniques : moving to convex combination");
coordinates.clear();
Octree_Delete(oct);
- fillNeumannBCS();
+ //fillNeumannBCS();
parametrize(ITERU,CONVEXCOMBINATION);
parametrize(ITERV,CONVEXCOMBINATION);
buildOct();
@@ -974,14 +974,14 @@ void GFaceCompound::parametrize(iterationStep step, typeOfMapping tom) const
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);
- }
- }
+// 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",
@@ -998,10 +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);
- }
+// 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);
@@ -1185,10 +1185,10 @@ void GFaceCompound::computeNormals() const
for(unsigned int i = 0; i < (*it)->triangles.size(); ++i){
MTriangle *t = (*it)->triangles[i];
t->getJacobian(0, 0, 0, J);
- // SVector3 n (J[2][0],J[2][1],J[2][2]);
- SVector3 d1(J[0][0], J[0][1], J[0][2]);
- SVector3 d2(J[1][0], J[1][1], J[1][2]);
- SVector3 n = crossprod(d1, d2);
+ SVector3 n (J[2][0],J[2][1],J[2][2]);
+ //SVector3 d1(J[0][0], J[0][1], J[0][2]);
+ //SVector3 d2(J[1][0], J[1][1], J[1][2]);
+ //SVector3 n = crossprod(d1, d2);
n.normalize();
for(int j = 0; j < 3; j++){
std::map<MVertex*, SVector3>::iterator itn = _normals.find(t->getVertex(j));
@@ -1253,7 +1253,7 @@ GPoint GFaceCompound::point(double par1, double par2) const
return lt->gf->point(pv.x(),pv.y());
}
- const bool LINEARMESH = false;
+ const bool LINEARMESH = true; //false
if(LINEARMESH){
@@ -1285,7 +1285,7 @@ GPoint GFaceCompound::point(double par1, double par2) const
w32 = dot(lt->v2 - lt->v3, n3);
w31 = dot(lt->v1 - lt->v3, n3);
w13 = dot(lt->v3 - lt->v1, n1);
- b210 = (2*lt->v1 + lt->v2 -w12*n1)*0.333;
+ b210 = (2*lt->v1 + lt->v2-w12*n1)*0.333;
b120 = (2*lt->v2 + lt->v1-w21*n2)*0.333;
b021 = (2*lt->v2 + lt->v3-w23*n2)*0.333;
b012 = (2*lt->v3 + lt->v2-w32*n3)*0.333;
@@ -1515,7 +1515,7 @@ bool GFaceCompound::checkAspectRatio() const
bool paramOK = true;
if(allNodes.empty()) buildAllNodes();
- double limit = 1.e15 ;
+ double limit = 1.e15;
double areaMax = 0.0;
int nb = 0;
std::list<GFace*>::const_iterator it = _compound.begin();
@@ -1548,7 +1548,8 @@ bool GFaceCompound::checkAspectRatio() const
SBoundingBox3d bboxH = bounds();
double H = norm(SVector3(bboxH.max(), bboxH.min()));
double D = getSizeBB(_U0);
- int split = (int)ceil(.7*H/D);
+ double eta = H/D;
+ int split = (int)ceil(eta/3);
nbSplit = std::max(split,2);
//printf("H=%g, D=%g split=%d nbSplit=%d \n", H, D, split, nbSplit);
Msg::Info("Partition geometry in N=%d parts", nbSplit);
diff --git a/Geo/discreteEdge.cpp b/Geo/discreteEdge.cpp
index 4100580492699e02033e0dfee575a302ade63358..01c2d9fd16ba287b907a4e45a7cc4ac7559e813a 100644
--- a/Geo/discreteEdge.cpp
+++ b/Geo/discreteEdge.cpp
@@ -400,7 +400,7 @@ GPoint discreteEdge::point(double par) const
MVertex *vB = lines[iEdge]->getVertex(0);
MVertex *vE = lines[iEdge]->getVertex(1);
- const bool LINEARMESH = false;
+ const bool LINEARMESH = true; //false;
if (LINEARMESH){
//linear Lagrange mesh
diff --git a/Mesh/multiscalePartition.cpp b/Mesh/multiscalePartition.cpp
index 84d282eaa39b48e4741123473d4cee1ccb41c00e..a4b7f8e6e3a48ec21fb429c280f1dd06838fb3b1 100644
--- a/Mesh/multiscalePartition.cpp
+++ b/Mesh/multiscalePartition.cpp
@@ -23,11 +23,9 @@ static void recur_connect (MVertex *v,
}
-static int connected_bounds (std::vector<MEdge> &edges)
+static int connected_bounds (std::vector<MEdge> &edges, std::vector<std::vector<MEdge> > &boundaries)
{
- 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++){
@@ -40,15 +38,17 @@ static int connected_bounds (std::vector<MEdge> &edges)
recur_connect (v2e.begin()->first,v2e,group,touched);
std::vector<MEdge> temp;
temp.insert(temp.begin(), group.begin(), group.end());
- regions.push_back(temp);
+ boundaries.push_back(temp);
for ( std::set<MVertex*>::iterator it = touched.begin() ; it != touched.end();++it)
v2e.erase(*it);
}
- return regions.size();
+ return boundaries.size();
}
-static int getGenus (std::vector<MElement *> &elements){
+
+static int getGenus (std::vector<MElement *> &elements,
+ std::vector<std::vector<MEdge> > &boundaries){
//We suppose MElements are simply connected
@@ -80,7 +80,7 @@ static int getGenus (std::vector<MElement *> &elements){
bEdges.erase(std::find(bEdges.begin(), bEdges.end(),me));
}
}
- nbBounds = connected_bounds(bEdges);
+ nbBounds = connected_bounds(bEdges, boundaries);
int genus = (int)(-poincare + 2 - nbBounds)/2;
//printf("************** partition has %d boundaries and genus =%d \n", nbBounds, genus);
@@ -89,31 +89,66 @@ static int getGenus (std::vector<MElement *> &elements){
}
-// static int getGeomAspectRatio (std::vector<MElement *> &elements){
+static int getAspectRatio (std::vector<MElement *> &elements,
+ std::vector<std::vector<MEdge> > &boundaries)
+{
+
+ std::set<MVertex*> vs;
+ for(unsigned int i = 0; i < elements.size(); i++){
+ MElement *e = elements[i];
+ for(int j = 0; j < e->getNumVertices(); j++){
+ vs.insert(e->getVertex(j));
+ }
+ }
+ SBoundingBox3d bb;
+ SOrientedBoundingBox obbox ;
+ std::vector<SPoint3> vertices;
+ for (std::set<MVertex* >::iterator it = vs.begin(); it != vs.end(); it++){
+ SPoint3 pt((*it)->x(),(*it)->y(), (*it)->z());
+ vertices.push_back(pt);
+ bb +=pt;
+ }
+
+// obb = SOrientedBoundingBox::buildOBB(vertices);
+// double H = obbox.getMaxSize();
+ double H = norm(SVector3(bb.max(), bb.min()));
+
+ double D = 0.0;
+ if (boundaries.size() == 0){
+ D=1.;
+ }
+ else{
+ for (int i=0; i< boundaries.size(); i++){
+ std::set<MVertex*> vb;
+ std::vector<MEdge> iBound = boundaries[i];
+ for (int j=0; j< iBound.size(); j++){
+ MEdge e = iBound[j];
+ vb.insert(e.getVertex(0));
+ vb.insert(e.getVertex(1));
+ }
+ std::vector<SPoint3> vBounds;
+ for (std::set<MVertex* >::iterator it = vb.begin(); it != vb.end(); it++){
+ SPoint3 pt((*it)->x(),(*it)->y(), (*it)->z());
+ vBounds.push_back(pt);
+ }
+ SOrientedBoundingBox obboxD = SOrientedBoundingBox::buildOBB(vBounds);
+ D = std::max(D, obboxD.getMaxSize());
+ }
+ }
+ int eta = (int)ceil(H/D);
+
+ return eta;
-// 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() );
+}
+static void getGenusAndRatio(std::vector<MElement *> &elements, int & genus, int &eta){
-// return eta;
+ std::vector<std::vector<MEdge> > boundaries;
+ genus = getGenus(elements, boundaries);
+ eta = getAspectRatio(elements, boundaries);
-// }
+ return;
+}
static void partitionRegions (std::vector<MElement*> &elements,
std::vector<std::vector<MElement*> > ®ions){
@@ -127,6 +162,7 @@ static void partitionRegions (std::vector<MElement*> &elements,
}
+
multiscalePartition::multiscalePartition (std::vector<MElement *> &elements,
meshPartitionOptions &_options){
@@ -144,6 +180,13 @@ multiscalePartition::multiscalePartition (std::vector<MElement *> &elements,
totalParts = assembleAllPartitions();
}
+void multiscalePartition:: setNumberOfPartitions(int &nbParts)
+{
+ options.num_partitions = nbParts;
+ if ( options.partitioner == 1){
+ options.mesh_dims[0] = nbParts;
+ }
+}
void multiscalePartition::partition(partitionLevel & level){
#if defined(HAVE_METIS) || defined(HAVE_CHACO)
@@ -161,20 +204,23 @@ void multiscalePartition::partition(partitionLevel & level){
nextLevel->region = i;
levels.push_back(nextLevel);
-
- int genus = getGenus(regions[i]);
+ int genus, eta;
+ getGenusAndRatio(regions[i], genus, eta);
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);
+ if (genus != 0 || eta > 4 ){
+ int nbParts = std::max(genus+2, (int)eta/4);
+ setNumberOfPartitions(nbParts);
+ Msg::Info("Multiscale partition, level %d region %d is %d-GENUS (ETA=%d) ---> partition %d parts",
+ nextLevel->recur,nextLevel->region, genus, eta, nbParts);
+
partition(*nextLevel);
}
else {
- Msg::Info("Multiscale Partition, level %d, region %d is ZERO-GENUS",
- nextLevel->recur,nextLevel->region);
+ Msg::Info("Multiscale Partition, level %d, region %d is ZERO-GENUS (ETA=%d)",
+ nextLevel->recur,nextLevel->region, eta);
}
}
@@ -185,7 +231,6 @@ void multiscalePartition::partition(partitionLevel & level){
int multiscalePartition::assembleAllPartitions(){
int nbParts = 1;
- //int nbLevel = options.num_partitions;
for (int i = 0; i< levels.size(); i++){
partitionLevel *iLevel = levels[i];
@@ -199,6 +244,13 @@ int multiscalePartition::assembleAllPartitions(){
}
}
+ Msg::Info("-----------------------------------------------------------");
+ Msg::Info("Multiscale Partition SUCCESSFULLY PERFORMED : %d parts", nbParts-1);
+ Msg::Info("-----------------------------------------------------------");
+
+
return nbParts-1;
+
+
}
diff --git a/Mesh/multiscalePartition.h b/Mesh/multiscalePartition.h
index 4d3db6d9e8e68d5d238c03958b1f81f21417b35d..0a5e1538d7fb1ecd7530a8ab0d58788020097d66 100644
--- a/Mesh/multiscalePartition.h
+++ b/Mesh/multiscalePartition.h
@@ -32,6 +32,7 @@ class multiscalePartition{
multiscalePartition(std::vector<MElement *> &elements,
meshPartitionOptions &options);
int assembleAllPartitions();
+ void setNumberOfPartitions(int &nbParts);
int getNumberOfParts(){return totalParts;}
};
diff --git a/contrib/Chaco/main/Gmsh_printf.cpp b/contrib/Chaco/main/Gmsh_printf.cpp
index 2f35dd581db66346c964c759a1c848ce1c484fa8..4ae2d6c75e348b41822f8491f2508c82bffe39d2 100644
--- a/contrib/Chaco/main/Gmsh_printf.cpp
+++ b/contrib/Chaco/main/Gmsh_printf.cpp
@@ -7,7 +7,7 @@
#include <cstring>
#include "GmshMessage.h"
-// Overload the printf statements in Chaco to write using Msg::Direct in gmsh
+// Overload the printf statements in Chaco to write using Msg::Debug in gmsh
extern "C" int Gmsh_printf(const char *fmt, ...)
{
@@ -23,17 +23,17 @@ extern "C" int Gmsh_printf(const char *fmt, ...)
char *p = std::strtok(str, "\n");
if(p) {
// If more than 1 leading '\n', print a blank line
- if(p - str > 1) Msg::Direct(" ");
+ if(p - str > 1) Msg::Debug(" ");
std::strcpy(buf, p);
- Msg::Direct(buf);
+ Msg::Debug(buf);
// New line for each interior '\n'
while((p = std::strtok(NULL, "\n"))) {
std::strcpy(buf, p);
- Msg::Direct(buf);
+ Msg::Debug(buf);
}
}
// If more than 1 trailing '\n', or only "\n" in the string, print a blank
// line.
- if(*last == '\n') Msg::Direct(" ");
+ if(*last == '\n') Msg::Debug(" ");
return 0;
}
diff --git a/contrib/Chaco/misc/countup_mesh.c b/contrib/Chaco/misc/countup_mesh.c
index 5d1ae7171ae11618b59833384a4e2f12e1723a8b..4e84eac5efce7ad435e0e78965369edfd9e088e6 100644
--- a/contrib/Chaco/misc/countup_mesh.c
+++ b/contrib/Chaco/misc/countup_mesh.c
@@ -241,22 +241,22 @@ int using_ewgts; /* are edge weights being used? */
}
ncuts /= 2;
nhops /= 2;
- printf("\n");
- printf(" Total Max/Set Min/Set\n");
- printf(" ----- ------- -------\n");
- printf("Set Size: %11d %11d %11d\n",
+ Gmsh_printf("\n");
+ Gmsh_printf(" Total Max/Set Min/Set\n");
+ Gmsh_printf(" ----- ------- -------\n");
+ Gmsh_printf("Set Size: %11d %11d %11d\n",
tot_size, max_size, min_size);
- printf("Edge Cuts: %11g %11g %11g\n",
+ Gmsh_printf("Edge Cuts: %11g %11g %11g\n",
ncuts, maxcuts, mincuts);
- printf("Mesh Hops: %11g %11g %11g\n",
+ Gmsh_printf("Mesh Hops: %11g %11g %11g\n",
nhops, maxhops, minhops);
- printf("Boundary Vertices: %11g %11g %11g\n",
+ Gmsh_printf("Boundary Vertices: %11g %11g %11g\n",
total_bdyvtxs, maxbdy, minbdy);
- printf("Boundary Vertex Hops: %11g %11g %11g\n",
+ Gmsh_printf("Boundary Vertex Hops: %11g %11g %11g\n",
bdyvtx_hops_tot, bdyvtx_hops_max, bdyvtx_hops_min);
- printf("Adjacent Sets: %11d %11d %11d\n",
+ Gmsh_printf("Adjacent Sets: %11d %11d %11d\n",
total_neighbors, maxneighbors, minneighbors);
- printf("Internal Vertices: %11d %11d %11d\n\n",
+ Gmsh_printf("Internal Vertices: %11d %11d %11d\n\n",
total_internal, max_internal, min_internal);
diff --git a/contrib/Chaco/submain/submain.c b/contrib/Chaco/submain/submain.c
index 24d9155c0008e9a9f68173b77c589f308f5c6973..ff34e0b65eb11b88aa659df239a2ccd026255fe2 100644
--- a/contrib/Chaco/submain/submain.c
+++ b/contrib/Chaco/submain/submain.c
@@ -172,7 +172,7 @@ long seed; /* for random graph mutations */
}
if (PRINT_HEADERS) {
- printf("\n\nStarting to partition ...\n\n");
+ Gmsh_printf("\n\nStarting to partition ...\n\n");
if (Output_File != NULL ) {
fprintf(Output_File,
"\n\nStarting to partition ... (residual, warning and error messages only)\n\n");