From a5880d08663440ff07a3e78b1b723d13843957d2 Mon Sep 17 00:00:00 2001
From: Emilie Marchandise <emilie.marchandise@uclouvain.be>
Date: Fri, 22 Jan 2010 16:46:30 +0000
Subject: [PATCH]
---
Geo/GFaceCompound.cpp | 53 ++++++++++----------
Geo/GModel.cpp | 62 +++--------------------
Geo/GModelIO_Mesh.cpp | 1 +
Geo/discreteEdge.cpp | 2 +-
Mesh/meshGFace.cpp | 5 +-
Mesh/meshPartition.cpp | 3 --
Mesh/multiscalePartition.cpp | 89 ++++++++++++++++++++++++++-------
Solver/linearSystemCSR.cpp | 2 +-
Solver/multiscaleLaplace.cpp | 96 ++++++++++++++++++++----------------
benchmarks/boolean/TORUS.geo | 15 +++++-
10 files changed, 180 insertions(+), 148 deletions(-)
diff --git a/Geo/GFaceCompound.cpp b/Geo/GFaceCompound.cpp
index 31d57fa222..8b9b4ab4e9 100644
--- a/Geo/GFaceCompound.cpp
+++ b/Geo/GFaceCompound.cpp
@@ -353,7 +353,7 @@ bool GFaceCompound::checkOrientation(int iter) const
int iterMax = 5;
if(!oriented && iter < iterMax){
if (iter == 0) Msg::Warning("*** Parametrization is NOT 1 to 1 : applying cavity checks.");
- Msg::Warning("*** Cavity Check - iter %d -",iter);
+ Msg::Debug("*** Cavity Check - iter %d -",iter);
one2OneMap();
return checkOrientation(iter+1);
}
@@ -451,7 +451,7 @@ bool GFaceCompound::parametrize() const
buildOct();
if (!checkOrientation(0)){
- Msg::Info("Parametrization failed using standard techniques : moving to convex combination");
+ Msg::Info("Parametrization switched to convex combination map");
coordinates.clear();
Octree_Delete(oct);
fillNeumannBCS();
@@ -719,7 +719,7 @@ GFaceCompound::GFaceCompound(GModel *m, int tag, std::list<GFace*> &compound,
if (!_lsys) {
#if defined(HAVE_PETSC)
_lsys = new linearSystemPETSc<double>;
-#elif defined(HAVE_TAUCS)
+#elif defined(_HAVE_TAUCS)
_lsys = new linearSystemCSRTaucs<double>;
#elif defined(HAVE_GMM)
linearSystemGmm<double> *_lsysb = new linearSystemGmm<double>;
@@ -1049,14 +1049,24 @@ void GFaceCompound::parametrize_conformal() const
}
MVertex *v1 = ordered[0];
- MVertex *v2 = ordered[1];
+ MVertex *v2 ;
+ double maxSize = 0.0;
+ for (int i=1; i< ordered.size(); i++){
+ MVertex *vi= ordered[i];
+ double dist = vi->distance(v1);
+ if (dist > maxSize){
+ v2 = vi;
+ maxSize = dist;
+ }
+ }
+
myAssembler.fixVertex(v1, 0, 1, 0);//0
myAssembler.fixVertex(v1, 0, 2, 0);//0
myAssembler.fixVertex(v2, 0, 1, 1);//1
myAssembler.fixVertex(v2, 0, 2, 0);//0
- //printf("Pinned vertex %g %g %g / %g %g %g \n", v1->x(), v1->y(), v1->z(), v2->x(), v2->y(), v2->z());
- //exit(1);
+ //printf("Pinned vertex %g %g %g / %g %g %g \n", v1->x(), v1->y(), v1->z(), v2->x(), v2->y(), v2->z());
+
std::list<GFace*>::const_iterator it = _compound.begin();
for( ; it != _compound.end(); ++it){
for(unsigned int i = 0; i < (*it)->triangles.size(); ++i){
@@ -1081,7 +1091,7 @@ void GFaceCompound::parametrize_conformal() const
simpleFunction<double> ONE(1.0);
- simpleFunction<double> MONE(1.0 );
+ simpleFunction<double> MONE(-1.0 );
laplaceTerm laplace1(model(), 1, &ONE);
laplaceTerm laplace2(model(), 2, &ONE);
crossConfTerm cross12(model(), 1, 2, &ONE);
@@ -1283,7 +1293,7 @@ GPoint GFaceCompound::point(double par1, double par2) const
// return lt->gf->point(pv.x(),pv.y());
// }
- const bool LINEARMESH = true; //false; //false
+ const bool LINEARMESH = true; //false
if(LINEARMESH){
@@ -1601,12 +1611,7 @@ bool GFaceCompound::checkTopology() const
double H = getSizeH();
double D = H;
- if (_interior_loops.size() > 0)
- D = getSizeBB(_U0);
-// for(std::list<std::list<GEdge*> >::const_iterator it = _interior_loops.begin(); it != _interior_loops.end(); it++){
-// double size = getSizeBB(*it);
-// D = std::min(D, size);
-// }
+ if (_interior_loops.size() > 0) D = getSizeBB(_U0);
int AR = (int) ceil(H/D);
if (G != 0 || Nb < 1){
@@ -1640,8 +1645,8 @@ bool GFaceCompound::checkAspectRatio() const
bool paramOK = true;
if(allNodes.empty()) buildAllNodes();
- double limit = 1.e12;
- double areaMax = 0.0;
+ double limit = 1.e-15;
+ double areaMin = 1.e15;
int nb = 0;
std::list<GFace*>::const_iterator it = _compound.begin();
for( ; it != _compound.end() ; ++it){
@@ -1662,22 +1667,19 @@ bool GFaceCompound::checkAspectRatio() const
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));
- if (1/a_2D > limit) nb++;
- areaMax = std::max(areaMax,1./a_2D);
- //printf("a3D/a2D=%g, AreaMax=%g\n", a_3D/a_2D, areaMax);
+ if (a_2D > limit) nb++;
+ areaMin = std::min(areaMin,a_2D);
}
}
- if (areaMax > limit && nb > 2) {
- Msg::Warning("Geometrical aspect ratio too high (1/area_2D=%g)", areaMax);
+ if (areaMin < limit && nb > 2) {
+ Msg::Warning("Geometrical aspect ratio too high (a_2D=%g)", areaMin);
SBoundingBox3d bboxH = bounds();
double H = getSizeH();
double D = getSizeBB(_U0);
double eta = H/D;
- int split = -2;
- //printf("H=%g, D=%g split=%d nbSplit=%d \n", H, D, split, nbSplit);
- Msg::Info("Partition geometry in N=%d parts", std::abs(nbSplit));
- paramOK = false;
+ int nbSplit = -2;
+ paramOK = true; //false;
}
else {
Msg::Debug("Geometrical aspect ratio is OK :-)");
@@ -1691,6 +1693,7 @@ bool GFaceCompound::checkAspectRatio() const
void GFaceCompound::coherenceNormals()
{
+ Msg::Info("Coherence Normals ");
std::map<MEdge, std::set<MTriangle*>, Less_Edge > edge2tris;
for(unsigned int i = 0; i < triangles.size(); i++){
MTriangle *t = triangles[i];
diff --git a/Geo/GModel.cpp b/Geo/GModel.cpp
index 0bd0b84f3d..dadb44e429 100644
--- a/Geo/GModel.cpp
+++ b/Geo/GModel.cpp
@@ -1078,6 +1078,9 @@ int GModel::removeDuplicateMeshVertices(double tolerance)
void GModel::createTopologyFromMesh()
{
+
+ Msg::Info("Creating topology from mesh");
+
// for each discreteRegion, create topology
std::vector<discreteRegion*> discRegions;
for(riter it = firstRegion(); it != lastRegion(); it++)
@@ -1085,7 +1088,7 @@ void GModel::createTopologyFromMesh()
discRegions.push_back((discreteRegion*) *it);
//EMI-FIX in case of createTopology for Volumes
- //all faces are set to the volume
+ //all faces are set to aeach volume
for (std::vector<discreteRegion*>::iterator it = discRegions.begin();
it != discRegions.end(); it++)
(*it)->setBoundFaces();
@@ -1096,63 +1099,14 @@ void GModel::createTopologyFromMesh()
if((*it)->geomType() == GEntity::DiscreteSurface)
discFaces.push_back((discreteFace*) *it);
-//EMI TODO
-//check for closed faces
-// for (std::vector<discreteFace*>::iterator itf = discFaces.begin(); itf != discFaces.end(); itf++){
-
-// std::list<MTriangles*> tris;
-// for (unsigned int i = 0; i < (*itf)->trinagles.size(); i++){
-// tris.push_back((*itf)->triangles[i]);
-// }
-
-// while (!tris.empty()) {
-// for (int i=0; i<3; i++) {
-// for (std::list<MTriangle*>::iterator it = tris.begin() ; it != segments.end(); ++it){
-// MEdge *e0 = (*it)->getEdge(0);
-// MEdge *e1 = (*it)->getEdge(1);
-// MEdge *e2 = (*it)->getEdge(2);
-// //printf("mline %d %d \n", v1->getNum(), v2->getNum());
-// std::list<MTriangle*>::iterator itp;
-// if ( v1 == vE ){
-// //printf("->push back this mline \n");
-// myLines.push_back(*it);
-// itp = it;
-// it++;
-// segments.erase(itp);
-// vE = v2;
-// i = -1;
-// }
-// else if ( v2 == vE){
-// //printf("->push back this mline \n");
-// 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;
-// //printf("not found VB=%d vE=%d\n", vB->getNum(), vE->getNum());
-// MVertex *temp = vB;
-// vB = vE;
-// vE = temp;
-// }
-// GFace *newGe = new discreteFace(GModel::current(), GModel::current()->maxFaceNum() + 1, 0, 0);
-// newGe->lines.insert(newGe->lines.end(), myLines.begin(), myLines.end());
-// GModel::current()->add(newGe);
-// }
-// }
-
-// }
+ //EMI TODO
+ //check for closed faces
- //create Topo From Faces
createTopologyFromFaces(discFaces);
- //create
+ //create old format (necessaray for boundary layers)
exportDiscreteGEOInternals();
+
}
void GModel::createTopologyFromFaces(std::vector<discreteFace*> &discFaces)
diff --git a/Geo/GModelIO_Mesh.cpp b/Geo/GModelIO_Mesh.cpp
index 0e623b4bd5..811e843869 100644
--- a/Geo/GModelIO_Mesh.cpp
+++ b/Geo/GModelIO_Mesh.cpp
@@ -839,6 +839,7 @@ int GModel::writeDistanceMSH(const std::string &name, double version, bool binar
//geometrical distance
+ //Compute distance for akk
std::map<MVertex*,double* > distance_map;
std::vector<SPoint3> pts;
std::vector<double> distances;
diff --git a/Geo/discreteEdge.cpp b/Geo/discreteEdge.cpp
index c5a2de40e8..e319a2a2b7 100644
--- a/Geo/discreteEdge.cpp
+++ b/Geo/discreteEdge.cpp
@@ -418,7 +418,7 @@ GPoint discreteEdge::point(double par) const
MVertex *vB = lines[iEdge]->getVertex(0);
MVertex *vE = lines[iEdge]->getVertex(1);
- const bool LINEARMESH = true; //false; //false;
+ const bool LINEARMESH = true; //false;
if (LINEARMESH){
//linear Lagrange mesh
diff --git a/Mesh/meshGFace.cpp b/Mesh/meshGFace.cpp
index 25f0f3cc86..e232b67521 100644
--- a/Mesh/meshGFace.cpp
+++ b/Mesh/meshGFace.cpp
@@ -1302,7 +1302,7 @@ bool checkMeshCompound(GFaceCompound *gf, std::list<GEdge*> &edges)
return isMeshed;
}
- bool correctParam = gf->parametrize();
+ bool correctParam = gf->parametrize();
if (!correctParam){
partitionAndRemesh((GFaceCompound*) gf);
@@ -1486,7 +1486,8 @@ void partitionAndRemesh(GFaceCompound *gf)
Msg::Info("*** Mesh of surface %d done by assembly remeshed faces", gf->tag());
Msg::Info("-----------------------------------------------------------");
- gf->coherenceNormals();
+
+ //gf->coherenceNormals();
gf->meshStatistics.status = GFace::DONE;
//CreateOutputFile("toto.msh", CTX::instance()->mesh.format);
diff --git a/Mesh/meshPartition.cpp b/Mesh/meshPartition.cpp
index 246ff3c795..d3ca6b7578 100644
--- a/Mesh/meshPartition.cpp
+++ b/Mesh/meshPartition.cpp
@@ -973,9 +973,6 @@ void assignPartitionBoundary(GModel *model,
ppe = new partitionEdge(model, -(int)pedges.size()-1, 0, 0, v2);
pedges.insert (ppe);
model->add(ppe);
- printf("*** Create partitionEdge %d (", ppe->tag());
- for (unsigned int i = 0; i < v2.size(); ++i) printf("%d ", v2[i]);
- printf(")\n");
}
else ppe = *it;
ppe->lines.push_back(new MLine (me.getVertex(0),me.getVertex(1)));
diff --git a/Mesh/multiscalePartition.cpp b/Mesh/multiscalePartition.cpp
index 3bdd227309..2df582db91 100644
--- a/Mesh/multiscalePartition.cpp
+++ b/Mesh/multiscalePartition.cpp
@@ -25,6 +25,25 @@ static void recur_connect(MVertex *v,
}
+
+// starting form a list of elements, returns
+// lists of lists that are all simply connected
+static void recur_connect_e (const MEdge &e,
+ std::multimap<MEdge,MElement*,Less_Edge> &e2e,
+ std::set<MElement*> &group,
+ std::set<MEdge,Less_Edge> &touched){
+ if (touched.find(e) != touched.end())return;
+ touched.insert(e);
+ for (std::multimap <MEdge,MElement*,Less_Edge>::iterator it = e2e.lower_bound(e);
+ it != e2e.upper_bound(e) ; ++it){
+ group.insert(it->second);
+ for (int i=0;i<it->second->getNumEdges();++i){
+ recur_connect_e (it->second->getEdge(i),e2e,group,touched);
+ }
+ }
+}
+
+
static int connected_bounds (std::vector<MEdge> &edges, std::vector<std::vector<MEdge> > &boundaries)
{
std::multimap<MVertex*,MEdge> v2e;
@@ -47,6 +66,27 @@ static int connected_bounds (std::vector<MEdge> &edges, std::vector<std::vector
return boundaries.size();
}
+//--------------------------------------------------------------
+static void connectedRegions (std::vector<MElement*> &elements,
+ std::vector<std::vector<MElement*> > ®ions)
+{
+ 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]));
+ }
+ }
+ while (!e2e.empty()){
+ std::set<MElement*> group;
+ std::set<MEdge,Less_Edge> touched;
+ recur_connect_e (e2e.begin()->first,e2e,group,touched);
+ std::vector<MElement*> temp;
+ temp.insert(temp.begin(), group.begin(), group.end());
+ regions.push_back(temp);
+ for ( std::set<MEdge>::iterator it = touched.begin() ; it != touched.end();++it)
+ e2e.erase(*it);
+ }
+}
static int getGenus (std::vector<MElement *> &elements,
std::vector<std::vector<MEdge> > &boundaries)
@@ -115,6 +155,7 @@ static int getAspectRatio(std::vector<MElement *> &elements,
double H = norm(SVector3(bb.max(), bb.min()));
double D = H;
+ if (boundaries.size() > 0.0 ) D = 0.0;
for (unsigned int i = 0; i < boundaries.size(); i++){
std::set<MVertex*> vb;
std::vector<MEdge> iBound = boundaries[i];
@@ -129,7 +170,7 @@ static int getAspectRatio(std::vector<MElement *> &elements,
vBounds.push_back(pt);
}
SOrientedBoundingBox obboxD = SOrientedBoundingBox::buildOBB(vBounds);
- D = std::min(D, obboxD.getMaxSize());
+ D = std::max(D, obboxD.getMaxSize());
}
int AR = (int)ceil(H/D);
@@ -146,11 +187,25 @@ static void getGenusAndRatio(std::vector<MElement *> &elements, int & genus, int
static void partitionRegions(std::vector<MElement*> &elements,
std::vector<std::vector<MElement*> > ®ions)
{
- for (unsigned int i = 0; i < elements.size(); ++i){
- MElement *e = elements[i];
- int part = e->getPartition();
- regions[part-1].push_back(e);
- }
+
+ for (unsigned int i = 0; i < elements.size(); ++i){
+ MElement *e = elements[i];
+ int part = e->getPartition();
+ regions[part-1].push_back(e);
+ }
+
+ std::vector<std::vector<MElement*> > allRegions;
+ for (unsigned int k = 0; k < regions.size(); ++k){
+ std::vector<std::vector<MElement*> > conRegions;
+ conRegions.clear();
+ connectedRegions (regions[k], conRegions);
+ for (int j=0; j< conRegions.size() ; j++)
+ allRegions.push_back(conRegions[j]);
+ }
+ regions.clear();
+ regions.resize(allRegions.size());
+ regions = allRegions;
+
}
static void printLevel(std::vector<MElement *> &elements, int recur, int region)
@@ -248,27 +303,27 @@ void multiscalePartition::partition(partitionLevel & level, int nbParts, typeOfP
int genus, AR;
getGenusAndRatio(regions[i], genus, AR);
+ printLevel (nextLevel->elements, nextLevel->recur,nextLevel->region);
+
if (genus < 0) {
Msg::Error("Genus partition is negative G=%d!", genus);
exit(1);
- }
-
- //printLevel (nextLevel->elements, nextLevel->recur,nextLevel->region);
-
+ }
+
if (genus != 0 ){
int nbParts = 2; //std::max(genus+2,2);
- Msg::Info("Multiscale part: level %d region %d is %d-GENUS (AR=%d) ---> MULTILEVEL partition %d parts",
+ 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);
}
else if (genus == 0 && AR > 3 ){
int nbParts = 2;
- Msg::Info("Multiscale part: level %d region %d is ZERO-GENUS (AR=%d) ---> LAPLACIAN partition %d parts",
+ Msg::Info("Mesh partition: level (%d-%d) is ZERO-GENUS (AR=%d) ---> LAPLACIAN partition %d parts",
nextLevel->recur,nextLevel->region, AR, nbParts);
partition(*nextLevel, nbParts, LAPLACIAN);
}
else {
- Msg::Info("*** Multiscale partition: level %d, region %d is ZERO-GENUS (AR=%d)",
+ Msg::Info("*** Mesh partition: level (%d-%d) is ZERO-GENUS (AR=%d)",
nextLevel->recur,nextLevel->region, AR);
}
@@ -279,7 +334,7 @@ void multiscalePartition::partition(partitionLevel & level, int nbParts, typeOfP
int multiscalePartition::assembleAllPartitions()
{
- int nbParts = 1;
+ int iPart = 1;
for (unsigned i = 0; i< levels.size(); i++){
partitionLevel *iLevel = levels[i];
@@ -287,11 +342,11 @@ int multiscalePartition::assembleAllPartitions()
for (unsigned j = 0; j < iLevel->elements.size(); j++){
MElement *e = iLevel->elements[j];
int part = e->getPartition();
- e->setPartition(nbParts);
+ e->setPartition(iPart);
}
- nbParts++;
+ iPart++;
}
}
- return nbParts - 1;
+ return iPart - 1;
}
diff --git a/Solver/linearSystemCSR.cpp b/Solver/linearSystemCSR.cpp
index 4bceac89ca..ec1739d006 100644
--- a/Solver/linearSystemCSR.cpp
+++ b/Solver/linearSystemCSR.cpp
@@ -352,7 +352,7 @@ int linearSystemCSRTaucs<double>::systemSolve()
options,
NULL);
double t2 = Cpu();
- Msg::Info("TAUCS has solved %d unknowns in %8.3f seconds", _b->size(), t2 - t1);
+ Msg::Debug("TAUCS has solved %d unknowns in %8.3f seconds", _b->size(), t2 - t1);
if (result != TAUCS_SUCCESS){
Msg::Error("Taucs Was Not Successfull %d", result);
}
diff --git a/Solver/multiscaleLaplace.cpp b/Solver/multiscaleLaplace.cpp
index 6d0f639489..b3ba1f1f7f 100644
--- a/Solver/multiscaleLaplace.cpp
+++ b/Solver/multiscaleLaplace.cpp
@@ -216,7 +216,8 @@ static void recur_compute_centers_ (double R, double a1, double a2,
SPoint2 PL (R*cos(a1),R*sin(a1));
SPoint2 PR (R*cos(a2),R*sin(a2));
centers.push_back(std::make_pair(PL,zero));
- centers.push_back(std::make_pair(PR,zero));
+ centers.push_back(std::make_pair(PR,zero));
+
for (int i=0;i<root->children.size();i++){
multiscaleLaplaceLevel* m = root->children[i];
centers.push_back(std::make_pair(m->center,m));
@@ -234,7 +235,7 @@ static void recur_compute_centers_ (double R, double a1, double a2,
connected_bounds(root->elements, boundaries);
//add the center of real holes ...
- if (root->children.size()==0 || boundaries.size()-1 != root->children.size()){
+ if (root->children.size()==0 ){//|| boundaries.size()-1 != root->children.size()){
for (int i = 0; i < boundaries.size(); i++){
std::vector<MEdge> me = boundaries[i];
SPoint2 c(0.0,0.0);
@@ -258,7 +259,6 @@ static void recur_compute_centers_ (double R, double a1, double a2,
//sort centers
std::sort(centers.begin(),centers.end());
- //printf("size centers =%d \n", centers.size());
for (int i=1;i<centers.size()-1;i++){
multiscaleLaplaceLevel* m1 = centers[i-1].second;
@@ -622,12 +622,32 @@ static void printLevel ( const char* fn,
}
//--------------------------------------------------------------
static double localSize(MElement *e, std::map<MVertex*,SPoint2> &solution){
+
SBoundingBox3d local;
for(unsigned int j = 0; j<e->getNumVertices(); ++j){
SPoint2 p = solution[e->getVertex(j)];
local += SPoint3(p.x(),p.y(),0.0);
}
- return local.max().distance(local.min());
+ return local.max().distance(local.min());
+
+// MVertex* v0 = e->getVertex(0);
+// MVertex* v1 = e->getVertex(1);
+// MVertex* v2 = e->getVertex(2);
+// double p0[3] = {v0->x(), v0->y(), v0->z()};
+// double p1[3] = {v1->x(), v1->y(), v1->z()};
+// double p2[3] = {v2->x(), v2->y(), v2->z()};
+// double a_3D = fabs(triangle_area(p0, p1, p2));
+// SPoint2 s1 = solution[v0];
+// SPoint2 s2 = solution[v1];
+// SPoint2 s3 = solution[v2];
+// double q0[3] = {s1.x(), s1.y(), 0.0};
+// double q1[3] = {s2.x(), s2.y(), 0.0};
+// double q2[3] = {s3.x(), s3.y(), 0.0};
+// double a_2D = fabs(triangle_area(q0, q1, q2));
+
+// return a_2D; //a_2D / a_3D;
+
+
}
//-------------------------------------------------------------
static void one2OneMap(std::vector<MElement *> &elements, std::map<MVertex*,SPoint2> &solution) {
@@ -700,8 +720,8 @@ static bool checkOrientation(std::vector<MElement *> &elements,
int iterMax = 1;
if(!oriented && iter < iterMax){
- if (iter == 0) Msg::Warning("*** Parametrization is NOT 1 to 1 : applying cavity checks.");
- Msg::Warning("*** Cavity Check - iter %d -",iter);
+ if (iter == 0) Msg::Debug("Parametrization is NOT 1 to 1 : applying cavity checks.");
+ Msg::Debug("Cavity Check - iter %d -",iter);
one2OneMap(elements, solution);
return checkOrientation(elements, solution, iter+1);
}
@@ -812,6 +832,7 @@ multiscaleLaplace::multiscaleLaplace (std::vector<MElement *> &elements, int iPa
void multiscaleLaplace::parametrize (multiscaleLaplaceLevel & level){
+ //Compute all nodes for the level
std::set<MVertex*> allNodes;
for(unsigned int i = 0; i < level.elements.size(); ++i){
MElement *e = level.elements[i];
@@ -824,10 +845,10 @@ void multiscaleLaplace::parametrize (multiscaleLaplaceLevel & level){
std::map<MVertex*,SPoint2> solution;
parametrize_method(level, allNodes, solution, HARMONIC);
if (!checkOrientation(level.elements, solution, 0)){
- Msg::Info("Parametrization failed using standard techniques : moving to convex combination");
+ Msg::Debug("Parametrization switched to convex combination");
parametrize_method(level, allNodes, solution, CONVEXCOMBINATION);
}
-
+
//Compute the bbox of the parametric space
SBoundingBox3d bbox;
for(std::set<MVertex *>::iterator itv = allNodes.begin(); itv !=allNodes.end() ; ++itv){
@@ -839,20 +860,16 @@ void multiscaleLaplace::parametrize (multiscaleLaplaceLevel & level){
//Check elements that are too small
std::vector<MElement*> tooSmall, goodSize;
-
for(unsigned int i = 0; i < level.elements.size(); ++i){
MElement *e = level.elements[i];
std::vector<SPoint2> localCoord;
double local_size = localSize(e,solution);
- if (local_size < 1.e-6 * global_size){
+ if (local_size < 1.e-6 * global_size)
tooSmall.push_back(e);
- }
- else{
- goodSize.push_back(e);
- }
+ else goodSize.push_back(e);
}
- //Only keep the right elements and nodes
+ //Only keep the connected elements vectors goodSize and tooSmall
std::vector<std::vector<MElement*> > regGoodSize;
connectedRegions (goodSize,regGoodSize);
int index=0;
@@ -868,30 +885,29 @@ void multiscaleLaplace::parametrize (multiscaleLaplaceLevel & level){
}
goodSize.clear();
for (int i=0;i< regGoodSize.size() ; i++){
- if (i == index)
- goodSize.insert(goodSize.begin(), regGoodSize[i].begin(), regGoodSize[i].end());
- else
- tooSmall.insert(tooSmall.begin(), regGoodSize[i].begin(), regGoodSize[i].end());
+ if (i == index) goodSize.insert(goodSize.begin(), regGoodSize[i].begin(), regGoodSize[i].end());
+ else tooSmall.insert(tooSmall.begin(), regGoodSize[i].begin(), regGoodSize[i].end());
}
+
level.elements = goodSize;
- std::vector<std::vector<MElement*> > regions, regions_;
+ //Add the not too small regions to the level.elements
+ std::vector<std::vector<MElement*> > regions_, regions ;
connectedRegions (tooSmall,regions_);
-
- //Remove small regions
for (int i=0;i< regions_.size() ; i++){
- bool region_has_really_small_elements = false;
+ bool really_small_elements = false;
for (int k=0; k<regions_[i].size() ; k++){
MElement *e = regions_[i][k];
double local_size = localSize(e,solution);
- if (local_size < 1.e-8 * global_size){
- region_has_really_small_elements = true;
- }
+ if (local_size < 1.e-8 * global_size)
+ really_small_elements = true;
}
- if(region_has_really_small_elements) regions.push_back(regions_[i]);
- else level.elements.insert(level.elements.begin(), regions_[i].begin(), regions_[i].end() );
+ if(really_small_elements && regions_[i].size() > 10) regions.push_back(regions_[i]);
+ else
+ level.elements.insert(level.elements.begin(), regions_[i].begin(), regions_[i].end() );
}
+ //Fill level.coordinates
std::set<MVertex*> goodSizev;
for(unsigned int i = 0; i < level.elements.size(); ++i){
MElement *e = level.elements[i];
@@ -900,23 +916,20 @@ void multiscaleLaplace::parametrize (multiscaleLaplaceLevel & level){
level.coordinates[e->getVertex(j)] = solution[e->getVertex(j)];
}
}
-
- // DEBUG
+
+ //Save multiscale meshes
char name[245];
sprintf(name,"multiscale_level%d_region%d_real.msh",level.recur, level.region);
printLevel (name,level.elements,0,2.0);
sprintf(name,"multiscale_level%d_region%d_param.msh",level.recur, level.region);
printLevel (name,level.elements,&level.coordinates,2.0);
- // END DEBUG
- if (regions.size() >0)
- Msg::Info("%d connected regions\n",regions.size());
+ //For every small region compute a new parametrization
+ Msg::Info("Level (%d-%d): %d connected small regions",level.recur,level.region, regions.size());
for (int i=0;i< regions.size() ; i++){
- // check nodes that appear both on too small and good size
- // and set it to the next level
- // maybe more than one level should be created here
std::set<MVertex*> tooSmallv;
+ tooSmallv.clear();
for (int k=0; k<regions[i].size() ; k++){
MElement *e = regions[i][k];
for(unsigned int j = 0; j<e->getNumVertices(); ++j){
@@ -927,7 +940,7 @@ void multiscaleLaplace::parametrize (multiscaleLaplaceLevel & level){
multiscaleLaplaceLevel *nextLevel = new multiscaleLaplaceLevel;
nextLevel->elements = regions[i];
nextLevel->recur = level.recur+1;
- nextLevel->region = i;
+ nextLevel->region = i;//Emi add something here
SBoundingBox3d smallB;
for(std::set<MVertex *>::iterator itv = tooSmallv.begin(); itv !=tooSmallv.end() ; ++itv){
SPoint2 p = solution[*itv];
@@ -944,16 +957,13 @@ void multiscaleLaplace::parametrize (multiscaleLaplaceLevel & level){
}
}
// recursively continue if tooSmall is not empty
- if (!tooSmall.empty()){
- Msg::Info("Multiscale Laplace, level %d region %d, %d too small",level.recur,level.region,tooSmall.size());
+ if (!tooSmallv.empty()){
+ Msg::Info("Level (%d-%d) Multiscale Laplace (reg[%d] = %d too small)",level.recur,level.region, i, tooSmallv.size());
level.children.push_back(nextLevel);
parametrize (*nextLevel);
}
- else {
- Msg::Info("Multiscale Laplace, level %d DONE",level.recur);
- delete nextLevel;
- }
}
+
}
void multiscaleLaplace::parametrize_method (multiscaleLaplaceLevel & level,
diff --git a/benchmarks/boolean/TORUS.geo b/benchmarks/boolean/TORUS.geo
index b890f1965e..e62b092c1c 100644
--- a/benchmarks/boolean/TORUS.geo
+++ b/benchmarks/boolean/TORUS.geo
@@ -1,2 +1,13 @@
-OCCShape("Torus",{0,0,0,0,0,1,10,9},"None");
-BoundingBox{0,0,0,10,10,10};
+Mesh.CharacteristicLengthFactor=0.1;
+
+X = 12;
+For I In {0:2}
+ OCCShape("Torus",{2*I*X,0,0,0,0,1,10,4},"Fuse");
+ OCCShape("Torus",{2*I*X,24,0,0,0,1,10,4},"Fuse");
+ OCCShape("Torus",{2*I*X,48,0,0,0,1,10,4},"Fuse");
+EndFor
+
+
+Compound Surface(100) = {1,2,3,4,5,6,7,8,9};
+
+
--
GitLab