diff --git a/Fltk/classificationEditor.cpp b/Fltk/classificationEditor.cpp
index 4a1249c7b3ec44e68728e7544b94c8bae00e80a7..d9ea4825c6beb185133665d5c6fa0f4692ae5c84 100644
--- a/Fltk/classificationEditor.cpp
+++ b/Fltk/classificationEditor.cpp
@@ -19,88 +19,6 @@
#include "discreteEdge.h"
#include "discreteFace.h"
-struct compareMLinePtr {
- bool operator () (MLine *l1, MLine *l2) const
- {
- static Less_Edge le;
- return le(l1->getEdge(0), l2->getEdge(0));
- }
-};
-
-static void recurClassify(MTri3 *t, GFace *gf,
- std::map<MLine*, GEdge*, compareMLinePtr> &lines,
- std::map<MTriangle*, GFace*> &reverse)
-{
- if(!t->isDeleted()){
- gf->triangles.push_back(t->tri());
- reverse[t->tri()] = gf;
- t->setDeleted(true);
- for(int i = 0; i < 3; i++){
- MTri3 *tn = t->getNeigh(i);
- if(tn){
- edgeXface exf(t, i);
- MLine ml(exf.v[0], exf.v[1]);
- std::map<MLine*, GEdge*, compareMLinePtr>::iterator it = lines.find(&ml);
- if(it == lines.end())
- recurClassify(tn, gf, lines, reverse);
- }
- }
- }
-}
-
-static GEdge *getNewModelEdge(GFace *gf1, GFace *gf2,
- std::map<std::pair<int, int>, GEdge*> &newEdges)
-{
- int t1 = gf1 ? gf1->tag() : -1;
- int t2 = gf2 ? gf2->tag() : -1;
- int i1 = std::min(t1, t2);
- int i2 = std::max(t1, t2);
-
- if(i1 == i2) return 0;
-
- std::map<std::pair<int, int>, GEdge*>::iterator it =
- newEdges.find(std::make_pair<int, int>(i1, i2));
- if(it == newEdges.end()){
- discreteEdge *ge = new discreteEdge
- (GModel::current(), GModel::current()->getMaxElementaryNumber(1) + 1, 0, 0);
- GModel::current()->add(ge);
- newEdges[std::make_pair<int, int>(i1, i2)] = ge;
- return ge;
- }
- else
- return it->second;
-}
-
-static void recurClassifyEdges(MTri3 *t,
- std::map<MTriangle*, GFace*> &reverse,
- std::map<MLine*, GEdge*, compareMLinePtr> &lines,
- std::set<MLine*> &touched,
- std::map<std::pair<int, int>, GEdge*> &newEdges)
-{
- if(!t->isDeleted()){
- t->setDeleted(true);
- GFace *gf1 = reverse[t->tri()];
- for(int i = 0; i < 3; i++){
- GFace *gf2 = 0;
- MTri3 *tn = t->getNeigh(i);
- if(tn)
- gf2 = reverse[tn->tri()];
- edgeXface exf(t, i);
- MLine ml(exf.v[0], exf.v[1]);
- std::map<MLine*, GEdge*, compareMLinePtr>::iterator it = lines.find(&ml);
- if(it != lines.end()){
- if(touched.find(it->first) == touched.end()){
- GEdge *ge = getNewModelEdge(gf1, gf2, newEdges);
- if(ge) ge->lines.push_back(it->first);
- touched.insert(it->first);
- }
- }
- if(tn)
- recurClassifyEdges(tn, reverse, lines, touched, newEdges);
- }
- }
-}
-
static void NoElementsSelectedMode(classificationEditor *e)
{
e->buttons[CLASS_BUTTON_SELECT_ELEMENTS]->activate();
@@ -372,118 +290,8 @@ static void select_surfaces_cb(Fl_Widget *w, void *data)
static void classify_cb(Fl_Widget *w, void *data)
{
classificationEditor *e = (classificationEditor*)data;
- std::map<MLine*, GEdge*, compareMLinePtr> lines;
- for(GModel::eiter it = GModel::current()->firstEdge();
- it != GModel::current()->lastEdge(); ++it){
- for(unsigned int i = 0; i < (*it)->lines.size();i++)
- lines[(*it)->lines[i]] = *it;
- }
-
- std::list<MTri3*> tris;
- {
- std::set<GFace*>::iterator it = e->faces.begin();
- while(it != e->faces.end()){
- GFace *gf = *it;
- for(unsigned int i = 0; i < gf->triangles.size(); i++)
- tris.push_back(new MTri3(gf->triangles[i], 0));
- gf->triangles.clear();
- ++it;
- }
- }
- if(tris.empty()) return;
-
- connectTriangles(tris);
-
- std::map<MTriangle*, GFace*> reverse;
- // color all triangles
- std::list<MTri3*> ::iterator it = tris.begin();
- while(it != tris.end()){
- if(!(*it)->isDeleted()){
- discreteFace *gf = new discreteFace
- (GModel::current(), GModel::current()->getMaxElementaryNumber(2) + 1);
- recurClassify(*it, gf, lines, reverse);
- GModel::current()->add(gf);
- }
- ++it;
- }
-
- // color some lines
- it = tris.begin();
- while(it != tris.end()){
- (*it)->setDeleted(false);
- ++it;
- }
-
- it = tris.begin();
-
- // classify edges that are bound by different GFaces
- std::map<std::pair<int, int>, GEdge*> newEdges;
- std::set<MLine*> touched;
- recurClassifyEdges(*it, reverse, lines, touched, newEdges);
-
- // check if new edges should not be splitted
- // splitted if composed of several open or closed edges
- for (std::map<std::pair<int, int>, GEdge*>::iterator it = newEdges.begin();
- it != newEdges.end() ; ++it){
- std::list<MLine*> segments;
- for(unsigned int i = 0; i < it->second->lines.size(); i++)
- segments.push_back(it->second->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;
- }
- GEdge *newGe = new discreteEdge
- (GModel::current(), GModel::current()->getMaxElementaryNumber(1) + 1, 0, 0);
- newGe->lines.insert(newGe->lines.end(), myLines.begin(), myLines.end());
- GModel::current()->add(newGe);
- }
- }
-
- for (std::map<std::pair<int, int>, GEdge*>::iterator it = newEdges.begin();
- it != newEdges.end(); ++it){
- GEdge *ge = it->second;
- GModel::current()->remove(ge);
- }
-
- while(it != tris.end()){
- delete *it;
- ++it;
- }
-
+ // classify faces with respect to coloured edges
+ GModel::current()->classifyFaces(e->faces);
// remove selected, but do not delete its elements
if(e->selected){
GModel::current()->remove(e->selected);
diff --git a/Fltk/statisticsWindow.cpp b/Fltk/statisticsWindow.cpp
index 546439453fe39fb445b54e0462a8b4a431a86ce2..576202d9ad0153764115415d27ec73efab7bc4e9 100644
--- a/Fltk/statisticsWindow.cpp
+++ b/Fltk/statisticsWindow.cpp
@@ -16,6 +16,7 @@
#include "Generator.h"
#include "Context.h"
#include "OS.h"
+#include "Field.h"
void statistics_cb(Fl_Widget *w, void *data)
{
@@ -204,15 +205,19 @@ void statisticsWindow::compute(bool elementQuality)
// meanAngle = meanAngle / count;
// printf("Angles = min=%g av=%g \n", minAngle, meanAngle);
//hack emi
-
+ /*
//Emi hack - MESH DEGREE VERTICES
std::vector<GEntity*> entities;
+ std::set<MEdge, Less_Edge> edges;
GModel::current()->getEntities(entities);
std::map<MVertex*, int > vert2Deg;
for(unsigned int i = 0; i < entities.size(); i++){
if(entities[i]->dim() < 2) continue;
for(unsigned int j = 0; j < entities[i]->getNumMeshElements(); j++){
MElement *e = entities[i]->getMeshElement(j);
+ for(unsigned int k = 0; k < e->getNumEdges(); k++){
+ edges.insert(e->getEdge(k));
+ }
for(unsigned int k = 0; k < e->getNumVertices(); k++){
MVertex *v = e->getVertex(k);
if (v->onWhat()->dim() < 2) continue;
@@ -243,6 +248,27 @@ void statisticsWindow::compute(bool elementQuality)
printf("Stats degree vertices: dMin=%d , dMax=%d, d4=%g \n", dMin, dMax, (double)d4/nbElems);
//end emi hack
+ FieldManager *fields = GModel::current()->getFields();
+ Field *f = fields->get(fields->background_field);
+ if(fields->background_field > 0){
+ std::set<MEdge, Less_Edge>::iterator it = edges.begin();
+ double sum = 0;
+ for (; it !=edges.end();++it){
+ MVertex *v0 = it->getVertex(0);
+ MVertex *v1 = it->getVertex(1);
+ double l = sqrt((v0->x()-v1->x())*(v0->x()-v1->x())+
+ (v0->y()-v1->y())*(v0->y()-v1->y())+
+ (v0->z()-v1->z())*(v0->z()-v1->z()));
+ double lf = (*f)(0.5*(v0->x()+v1->x()), 0.5*(v0->y()+v1->y()),0.5*(v0->z()+v1->z()),v0->onWhat());
+ double e = (l>lf) ? lf/l : l/lf;
+ sum += e - 1.0;
+ }
+ double tau = exp ((1./edges.size()) * sum);
+ printf("nedegs = %d tau = %g\n",edges.size(),tau);
+ }
+ */
+
+
int num = 0;
static double s[50];
static char label[50][256];
diff --git a/Geo/GModel.cpp b/Geo/GModel.cpp
index a6c8b05cb93d7fa789c5f33c8eba92616f255b88..4197747934f6fc2eb03e126659200f241e523da1 100644
--- a/Geo/GModel.cpp
+++ b/Geo/GModel.cpp
@@ -35,6 +35,7 @@
#include "MVertexPositionSet.h"
#include "OpenFile.h"
#include "CreateFile.h"
+#include "meshGFaceOptimize.h"
#if defined(HAVE_MESH)
#include "Field.h"
@@ -1993,6 +1994,339 @@ void GModel::insertRegion(GRegion *r)
regions.insert(r);
}
+// given a set of mesh edges, build GFaces that separate those
+// edges.
+
+struct compareMLinePtr {
+ bool operator () (MLine *l1, MLine *l2) const
+ {
+ static Less_Edge le;
+ return le(l1->getEdge(0), l2->getEdge(0));
+ }
+};
+
+
+static GEdge *getNewModelEdge(GFace *gf1, GFace *gf2,
+ std::map<std::pair<int, int>, GEdge*> &newEdges)
+{
+ int t1 = gf1 ? gf1->tag() : -1;
+ int t2 = gf2 ? gf2->tag() : -1;
+ int i1 = std::min(t1, t2);
+ int i2 = std::max(t1, t2);
+
+ if(i1 == i2) return 0;
+
+ std::map<std::pair<int, int>, GEdge*>::iterator it =
+ newEdges.find(std::make_pair<int, int>(i1, i2));
+ if(it == newEdges.end()){
+ discreteEdge *ge = new discreteEdge
+ (GModel::current(), GModel::current()->getMaxElementaryNumber(1) + 1, 0, 0);
+ GModel::current()->add(ge);
+ newEdges[std::make_pair<int, int>(i1, i2)] = ge;
+ return ge;
+ }
+ else
+ return it->second;
+}
+
+static void recurClassifyEdges(MTri3 *t,
+ std::map<MTriangle*, GFace*> &reverse,
+ std::map<MLine*, GEdge*, compareMLinePtr> &lines,
+ std::set<MLine*> &touched,
+ std::set<MTri3*> &trisTouched,
+ std::map<std::pair<int, int>, GEdge*> &newEdges)
+{
+ if(!t->isDeleted()){
+ trisTouched.erase(t);
+ t->setDeleted(true);
+ GFace *gf1 = reverse[t->tri()];
+ for(int i = 0; i < 3; i++){
+ GFace *gf2 = 0;
+ MTri3 *tn = t->getNeigh(i);
+ if(tn)
+ gf2 = reverse[tn->tri()];
+ edgeXface exf(t, i);
+ MLine ml(exf.v[0], exf.v[1]);
+ std::map<MLine*, GEdge*, compareMLinePtr>::iterator it = lines.find(&ml);
+ if(it != lines.end()){
+ if(touched.find(it->first) == touched.end()){
+ GEdge *ge = getNewModelEdge(gf1, gf2, newEdges);
+ if(ge) ge->lines.push_back(it->first);
+ touched.insert(it->first);
+ }
+ }
+ if(tn)
+ recurClassifyEdges(tn, reverse, lines, touched, trisTouched,newEdges);
+ }
+ }
+}
+
+static void recurClassify(MTri3 *t, GFace *gf,
+ std::map<MLine*, GEdge*, compareMLinePtr> &lines,
+ std::map<MTriangle*, GFace*> &reverse)
+{
+ if(!t->isDeleted()){
+ gf->triangles.push_back(t->tri());
+ reverse[t->tri()] = gf;
+ t->setDeleted(true);
+ for(int i = 0; i < 3; i++){
+ MTri3 *tn = t->getNeigh(i);
+ if(tn){
+ edgeXface exf(t, i);
+ MLine ml(exf.v[0], exf.v[1]);
+ std::map<MLine*, GEdge*, compareMLinePtr>::iterator it = lines.find(&ml);
+ if(it == lines.end())
+ recurClassify(tn, gf, lines, reverse);
+ }
+ }
+ }
+}
+
+void GModel::detectEdges(double _tresholdAngle) {
+ e2t_cont adj;
+ std::vector<MTriangle*> elements;
+ std::vector<edge_angle> edges_detected, edges_lonly;
+ for(GModel::fiter it = GModel::current()->firstFace();
+ it != GModel::current()->lastFace(); ++it)
+ elements.insert(elements.end(), (*it)->triangles.begin(),
+ (*it)->triangles.end());
+ buildEdgeToTriangle(elements, adj);
+ buildListOfEdgeAngle(adj, edges_detected, edges_lonly);
+ GEdge *selected = new discreteEdge
+ (this, getMaxElementaryNumber(1) + 1, 0, 0);
+ add(selected);
+
+ for(unsigned int i = 0; i < edges_detected.size(); i++){
+ edge_angle ea = edges_detected[i];
+ if(ea.angle <= _tresholdAngle) break;
+ selected->lines.push_back(new MLine(ea.v1, ea.v2));
+ }
+
+ for(unsigned int i = 0 ; i < edges_lonly.size(); i++){
+ edge_angle ea = edges_lonly[i];
+ selected->lines.push_back(new MLine(ea.v1, ea.v2));
+ }
+ std::set<GFace*> _temp;
+ _temp.insert(faces.begin(),faces.end());
+ classifyFaces(_temp);
+ remove(selected);
+ // delete selected;
+}
+
+void GModel::classifyFaces(std::set<GFace*> &_faces) {
+
+ std::map<MLine*, GEdge*, compareMLinePtr> lines;
+
+ for(GModel::eiter it = GModel::current()->firstEdge();
+ it != GModel::current()->lastEdge(); ++it){
+ for(unsigned int i = 0; i < (*it)->lines.size();i++)
+ lines[(*it)->lines[i]] = *it;
+ }
+
+ std::map<MTriangle*, GFace*> reverse_old;
+ std::list<MTri3*> tris;
+ {
+ std::set<GFace*>::iterator it = _faces.begin();
+ while(it != _faces.end()){
+ GFace *gf = *it;
+ for(unsigned int i = 0; i < gf->triangles.size(); i++){
+ tris.push_back(new MTri3(gf->triangles[i], 0));
+ reverse_old[gf->triangles[i]] = gf;
+ }
+ gf->triangles.clear();
+ gf->mesh_vertices.clear();
+ ++it;
+ }
+ }
+ if(tris.empty()) return;
+
+ connectTriangles(tris);
+
+ std::map<MTriangle*, GFace*> reverse;
+ std::multimap<GFace*, GFace*> replacedBy;
+ // color all triangles
+ std::list<MTri3*> ::iterator it = tris.begin();
+ std::list<GFace*> newf;
+ while(it != tris.end()){
+ if(!(*it)->isDeleted()){
+ discreteFace *gf = new discreteFace
+ (GModel::current(), GModel::current()->getMaxElementaryNumber(2) + 1);
+ recurClassify(*it, gf, lines, reverse);
+ GModel::current()->add(gf);
+ newf.push_back(gf);
+
+ for (int i=0;i<gf->triangles.size();i++){
+ replacedBy.insert(std::make_pair(reverse_old[gf->triangles[i]],gf));
+ }
+ }
+ ++it;
+ }
+
+ // now we have all faces coloured. If some regions were existing, replace
+ // their faces by the new ones
+
+ for (riter rit = firstRegion(); rit != lastRegion(); ++rit){
+ std::list<GFace *> _xfaces = (*rit)->faces();
+ std::set<GFace *> _newFaces;
+ for (std::list<GFace *>::iterator itf = _xfaces.begin(); itf != _xfaces.end(); ++itf){
+ std::multimap<GFace*, GFace*>::iterator itLow = replacedBy.lower_bound(*itf);
+ std::multimap<GFace*, GFace*>::iterator itUp = replacedBy.upper_bound(*itf);
+ for (; itLow != itUp; ++itLow)
+ _newFaces.insert(itLow->second);
+ }
+ std::list<GFace *> _temp;
+ _temp.insert(_temp.begin(),_newFaces.begin(),_newFaces.end());
+ (*rit)->set(_temp);
+ }
+
+ // color some lines
+ it = tris.begin();
+ while(it != tris.end()){
+ (*it)->setDeleted(false);
+ ++it;
+ }
+
+ // classify edges that are bound by different GFaces
+ std::map<std::pair<int, int>, GEdge*> newEdges;
+ std::set<MLine*> touched;
+ std::set<MTri3*> trisTouched;
+ // bug fix : multiply connected domains
+
+ trisTouched.insert(tris.begin(),tris.end());
+ while(!trisTouched.empty())
+ recurClassifyEdges(*trisTouched.begin(), reverse, lines, touched, trisTouched,newEdges);
+
+ // printf("%d new edges\n",newEdges.size());
+
+ std::map<discreteFace*,std::vector<int> > newFaceTopology;
+
+ // check if new edges should not be splitted
+ // splitted if composed of several open or closed edges
+
+ std::map<MVertex*,GVertex*> modelVertices;
+
+ for (std::map<std::pair<int, int>, GEdge*>::iterator ite = newEdges.begin();
+ ite != newEdges.end() ; ++ite){
+ std::list<MLine*> allSegments;
+ for(unsigned int i = 0; i < ite->second->lines.size(); i++)
+ allSegments.push_back(ite->second->lines[i]);
+
+ while (!allSegments.empty()) {
+ std::list<MLine*> segmentsForThisDiscreteEdge;
+ MVertex *vB = (*allSegments.begin())->getVertex(0);
+ MVertex *vE = (*allSegments.begin())->getVertex(1);
+ segmentsForThisDiscreteEdge.push_back(*allSegments.begin());
+ allSegments.erase(allSegments.begin());
+ while(1){
+ bool found = false;
+ for (std::list<MLine*>::iterator it = allSegments.begin();it != allSegments.end(); ++it){
+ MVertex *v1 = (*it)->getVertex(0);
+ MVertex *v2 = (*it)->getVertex(1);
+ if (v1 == vE || v2 == vE){
+ segmentsForThisDiscreteEdge.push_back(*it);
+ if (v2 == vE)(*it)->revert();
+ vE = (v1 == vE) ? v2 : v1;
+ found = true;
+ allSegments.erase(it);
+ break;
+ }
+ if (v1 == vB || v2 == vB){
+ segmentsForThisDiscreteEdge.push_front(*it);
+ if (v1 == vB)(*it)->revert();
+ vB = (v1 == vB) ? v2 : v1;
+ found = true;
+ allSegments.erase(it);
+ break;
+ }
+ }
+ if (vE == vB)break;
+ if (!found)break;
+ }
+
+ std::map<MVertex*,GVertex*>::iterator itMV = modelVertices.find(vB);
+ if (itMV == modelVertices.end()){
+ GVertex *newGv = new discreteVertex(GModel::current(),
+ GModel::current()->getMaxElementaryNumber(0) + 1);
+ newGv->mesh_vertices.push_back(vB);
+ vB->setEntity(newGv);
+ newGv->points.push_back(new MPoint(vB));
+ GModel::current()->add(newGv);
+ modelVertices[vB] = newGv;
+ }
+ itMV = modelVertices.find(vE);
+ if (itMV == modelVertices.end()){
+ GVertex *newGv = new discreteVertex(GModel::current(),
+ GModel::current()->getMaxElementaryNumber(0) + 1);
+ newGv->mesh_vertices.push_back(vE);
+ newGv->points.push_back(new MPoint(vE));
+ vE->setEntity(newGv);
+ GModel::current()->add(newGv);
+ modelVertices[vE] = newGv;
+ }
+
+ GEdge *newGe = new discreteEdge
+ (GModel::current(), GModel::current()->getMaxElementaryNumber(1) + 1, modelVertices[vB], modelVertices[vE]);
+ newGe->lines.insert(newGe->lines.end(), segmentsForThisDiscreteEdge.begin(), segmentsForThisDiscreteEdge.end());
+
+ // printf("discrete edge %d\n",newGe->tag());
+ for (std::list<MLine*>::iterator itL = segmentsForThisDiscreteEdge.begin();
+ itL != segmentsForThisDiscreteEdge.end(); ++itL){
+ if((*itL)->getVertex(0)->onWhat()->dim() != 0){
+ newGe->mesh_vertices.push_back((*itL)->getVertex(0));
+ (*itL)->getVertex(0)->setEntity(newGe);
+ }
+ // printf("%d %d \n",(*itL)->getVertex(0)->getNum(),(*itL)->getVertex(1)->getNum());
+ }
+
+ GModel::current()->add(newGe);
+ discreteFace *gf1 = dynamic_cast<discreteFace*> (GModel::current()->getFaceByTag(ite->first.first));
+ discreteFace *gf2 = dynamic_cast<discreteFace*> (GModel::current()->getFaceByTag(ite->first.second));
+ if (gf1)newFaceTopology[gf1].push_back(newGe->tag());
+ if (gf2)newFaceTopology[gf2].push_back(newGe->tag());
+ }
+ }
+
+ std::map<discreteFace*,std::vector<int> >::iterator itFT = newFaceTopology.begin();
+ for (;itFT != newFaceTopology.end();++itFT){
+ itFT->first->setBoundEdges(itFT->second);
+ }
+
+
+ for (std::map<std::pair<int, int>, GEdge*>::iterator it = newEdges.begin();
+ it != newEdges.end(); ++it){
+ GEdge *ge = it->second;
+ GModel::current()->remove(ge);
+ // delete ge;
+ }
+
+ it = tris.begin();
+ while(it != tris.end()){
+ delete *it;
+ ++it;
+ }
+
+ // delete empty mesh faces and reclasssify
+ std::set<GFace*, GEntityLessThan> fac = faces;
+ for (fiter fit = fac.begin() ; fit !=fac.end() ; ++fit){
+ std::set<MVertex *> _verts;
+ (*fit)->mesh_vertices.clear();
+ for (int i=0;i<(*fit)->triangles.size();i++){
+ for (int j=0;j<3;j++){
+ if ((*fit)->triangles[i]->getVertex(j)->onWhat()->dim() > 1){
+ (*fit)->triangles[i]->getVertex(j)->setEntity(*fit);
+ _verts.insert((*fit)->triangles[i]->getVertex(j));
+ }
+ }
+ }
+ if ((*fit)->triangles.size())
+ (*fit)->mesh_vertices.insert((*fit)->mesh_vertices.begin(),_verts.begin(),_verts.end());
+ else
+ remove(*fit);
+ }
+}
+
+
+
#include "Bindings.h"
void GModel::registerBindings(binding *b)
@@ -2155,6 +2489,10 @@ void GModel::registerBindings(binding *b)
cm->setDescription("build the topology for a given mesh.");
cm->setArgNames(NULL);
+ cm = cb->addMethod("detectEdges", &GModel::detectEdges);
+ cm->setDescription(" use an angle threshold to tag edges.");
+ cm->setArgNames("angle",NULL);
+
cm = cb->addMethod("createBoundaryLayer", &GModel::createBoundaryLayer);
cm->setDescription("create a boundary layer using a given field for the "
"extrusion height.");
diff --git a/Geo/GModel.h b/Geo/GModel.h
index 8a1297ea358f7e0bda994fcde11720086cfb54d0..1b491bfd78adb10da6c3b04b8b66b9a723f57d48 100644
--- a/Geo/GModel.h
+++ b/Geo/GModel.h
@@ -362,6 +362,11 @@ class GModel
// mesh the model
int mesh(int dimension);
+ // reclassify a mesh i.e. use an angle threshold to tag edges faces and regions
+ void detectEdges(double _tresholdAngle);
+ void classifyFaces(std::set<GFace*> &_faces);
+
+
// glue entities in the model (assume a tolerance eps and merge
// vertices that are too close, then merge edges, faces and
// regions). Warning: the gluer changes the geometric model, so that
diff --git a/Geo/GModelIO_Mesh.cpp b/Geo/GModelIO_Mesh.cpp
index e6d1f18a8b6fb37a644c17dae326b315a8cbde96..1477f38978b16d8ba0a13b538886f42700714c4a 100644
--- a/Geo/GModelIO_Mesh.cpp
+++ b/Geo/GModelIO_Mesh.cpp
@@ -1077,9 +1077,14 @@ int GModel::readSTL(const std::string &name, double tolerance)
}
Msg::Info("%d facets in solid %d", points[i].size() / 3, i);
// create face
- GFace *face = new discreteFace(this, getMaxElementaryNumber(2) + 1);
+ GFace *face = new discreteFace(this, getMaxElementaryNumber(2) + 1);
+ GRegion *region = new GRegion(this, getMaxElementaryNumber(3) + 1);
faces.push_back(face);
add(face);
+ add(region);
+ std::list<GFace*> _temp;
+ _temp.push_back(face);
+ region->set(_temp);
}
// create triangles using unique vertices
@@ -1091,6 +1096,8 @@ int GModel::readSTL(const std::string &name, double tolerance)
points[i][j].z()));
MVertexPositionSet pos(vertices);
+ std::set<MFace,Less_Face> unique;
+ int nbDuplic = 0;
for(unsigned int i = 0; i < points.size(); i ++){
for(unsigned int j = 0; j < points[i].size(); j += 3){
MVertex *v[3];
@@ -1100,9 +1107,17 @@ int GModel::readSTL(const std::string &name, double tolerance)
double z = points[i][j + k].z();
v[k] = pos.find(x, y, z, eps);
}
- faces[i]->triangles.push_back(new MTriangle(v[0], v[1], v[2]));
+ MFace mf (v[0],v[1],v[2]);
+ if (unique.find(mf) == unique.end()){
+ faces[i]->triangles.push_back(new MTriangle(v[0], v[1], v[2]));
+ unique.insert(mf);
+ }
+ else{
+ nbDuplic++;
+ }
}
}
+ if (nbDuplic) Msg::Warning("%d Duplicate triangle in STL file read",nbDuplic);
_associateEntityWithMeshVertices();
_storeVerticesInEntities(vertices); // will delete unused vertices
diff --git a/Geo/GModelIO_OCC.cpp b/Geo/GModelIO_OCC.cpp
index 8e1d5eb302ade22fe2f925cbf3e0ec0096843da2..cc14fe9075a5159cf0e3791f88dca14511db0e55 100644
--- a/Geo/GModelIO_OCC.cpp
+++ b/Geo/GModelIO_OCC.cpp
@@ -501,6 +501,7 @@ GRegion* OCC_Internals::addRegionToModel(GModel *model, TopoDS_Solid region){
void OCC_Internals::buildGModel(GModel *model)
{
+
// building geom vertices
int nvertices = vmap.Extent();
for(int i = 1; i <= nvertices; i++){
@@ -508,6 +509,7 @@ void OCC_Internals::buildGModel(GModel *model)
if (!getOCCVertexByNativePtr(model,TopoDS::Vertex(vmap(i))))
model->add(new OCCVertex(model,num , TopoDS::Vertex(vmap(i))));
}
+
// building geom edges
int nedges = emap.Extent();
for(int i = 1; i <= nedges; i++){
@@ -517,6 +519,7 @@ void OCC_Internals::buildGModel(GModel *model)
if (!getOCCEdgeByNativePtr(model,TopoDS::Edge(emap(i))))
model->add(new OCCEdge(model, TopoDS::Edge(emap(i)), num, model->getVertexByTag(i1), model->getVertexByTag(i2)));
}
+
// building geom faces
int nfaces = fmap.Extent();
for(int i = 1; i <= nfaces; i++){
diff --git a/Geo/GRegion.cpp b/Geo/GRegion.cpp
index 3741d99f2aabe20a8021da33942e45b6ec5d74f8..fa43377ae45a9a9e73005682f3c9b30954bdd9f2 100644
--- a/Geo/GRegion.cpp
+++ b/Geo/GRegion.cpp
@@ -299,6 +299,90 @@ void GRegion::replaceFaces (std::list<GFace*> &new_faces)
l_dirs = newdirs;
}
+double GRegion::computeSolidProperties (std::vector<double> cg,
+ std::vector<double> inertia){
+ std::list<GFace*>::iterator it = l_faces.begin();
+ std::list<int>::iterator itdir = l_dirs.begin();
+ double volumex = 0;
+ double volumey = 0;
+ double volumez = 0;
+ double surface = 0;
+ cg[0] = cg[1] = cg[2] = 0.0;
+ for ( ; it != l_faces.end(); ++it,++itdir){
+ printf("face %d dir %d %d elements\n",(*it)->tag(),*itdir,(*it)->triangles.size());
+ for (int i=0;i<(*it)->triangles.size();++i){
+ MTriangle *e = (*it)->triangles[i];
+ // MElement *e = (*it)->getMeshElement(i);
+ int npt;
+ IntPt *pts;
+ e->getIntegrationPoints (2*(e->getPolynomialOrder()-1)+3, &npt, &pts);
+ for (int j=0;j<npt;j++){
+ SPoint3 pt;
+ // compute x,y,z of the integration point
+ e->pnt(pts[j].pt[0],pts[j].pt[1],pts[j].pt[2],pt);
+ double jac[3][3];
+ // compute normal
+ double detJ = e->getJacobian(pts[j].pt[0],pts[j].pt[1],pts[j].pt[2],jac);
+ SVector3 n (jac[2][0],jac[2][1],jac[2][2]);
+ n.normalize();
+ n *= (double)*itdir;
+ surface += detJ*pts[j].weight;
+ volumex += detJ*n.x()*pt.x()*pts[j].weight;
+ volumey += detJ*n.y()*pt.y()*pts[j].weight;
+ volumez += detJ*n.z()*pt.z()*pts[j].weight;
+ cg[0] += detJ*n.x()*(pt.x()*pt.x())*pts[j].weight*0.5;
+ cg[1] += detJ*n.y()*(pt.y()*pt.y())*pts[j].weight*0.5;
+ cg[2] += detJ*n.z()*(pt.z()*pt.z())*pts[j].weight*0.5;
+ }
+ }
+ }
+
+ printf("%g -- %g %g %g\n",surface,volumex,volumey,volumez);
+
+ double volume = volumex;
+
+ cg[0]/=volume;
+ cg[1]/=volume;
+ cg[2]/=volume;
+
+ it = l_faces.begin();
+ itdir = l_dirs.begin();
+ inertia[0] =
+ inertia[1] =
+ inertia[2] =
+ inertia[3] =
+ inertia[4] =
+ inertia[5] = 0.0;
+
+ for ( ; it != l_faces.end(); ++it,++itdir){
+ for (int i=0;i<(*it)->getNumMeshElements();++i){
+ MElement *e = (*it)->getMeshElement(i);
+ int npt;
+ IntPt *pts;
+ e->getIntegrationPoints (2*(e->getPolynomialOrder()-1)+3, &npt, &pts);
+ for (int j=0;j<npt;j++){
+ SPoint3 pt;
+ // compute x,y,z of the integration point
+ e->pnt(pts[j].pt[0],pts[j].pt[1],pts[j].pt[2],pt);
+ double jac[3][3];
+ // compute normal
+ double detJ = e->getJacobian(pts[j].pt[0],pts[j].pt[1],pts[j].pt[2],jac);
+ SVector3 n (jac[2][0],jac[2][1],jac[2][2]);
+ n *= (double)*itdir;
+ inertia[0] += pts[j].weight*detJ*n.x()*(pt.x()-cg[0])*(pt.x()-cg[0])*(pt.x()-cg[0])/3.0;
+ inertia[1] += pts[j].weight*detJ*n.y()*(pt.y()-cg[1])*(pt.y()-cg[1])*(pt.y()-cg[1])/3.0;
+ inertia[2] += pts[j].weight*detJ*n.z()*(pt.z()-cg[2])*(pt.z()-cg[2])*(pt.z()-cg[2])/3.0;
+ inertia[3] += pts[j].weight*detJ*n.x()*(pt.y()-cg[1])*(pt.x()-cg[0])*(pt.x()-cg[0])/3.0;
+ inertia[4] += pts[j].weight*detJ*n.x()*(pt.z()-cg[2])*(pt.x()-cg[0])*(pt.x()-cg[0])/3.0;
+ inertia[5] += pts[j].weight*detJ*n.y()*(pt.z()-cg[2])*(pt.y()-cg[1])*(pt.y()-cg[1])/3.0;
+ }
+ }
+ }
+ return volume;
+}
+
+
+
void GRegion::addPrism(MPrism *p) {
prisms.push_back(p);
}
@@ -317,4 +401,7 @@ void GRegion::registerBindings(binding *b)
cm = cb->addMethod("addPrism", &GRegion::addPrism);
cm->setDescription("insert a prism mesh element");
cm->setArgNames("prism", NULL);
+ cm = cb->addMethod("computeSolidProperties", &GRegion::computeSolidProperties);
+ cm->setDescription("returns the volume and computes the center of gravity and tensor of inertia of the volume (requires a surface mesh)");
+ cm->setArgNames("cg","inertia", NULL);
}
diff --git a/Geo/GRegion.h b/Geo/GRegion.h
index 705361afee2a9989c5c5614993b4374f9436fc58..d00fe5ddc877a9e3bec56948460ed7dd29c3a593 100644
--- a/Geo/GRegion.h
+++ b/Geo/GRegion.h
@@ -107,6 +107,10 @@ class GRegion : public GEntity {
// replace edges (gor gluing)
void replaceFaces (std::list<GFace*> &);
+ // compute volume, moment of intertia and center of gravity
+ double computeSolidProperties (std::vector<double> cg,
+ std::vector<double> inertia);
+
static void registerBindings(binding *b);
};
diff --git a/Geo/Geo.cpp b/Geo/Geo.cpp
index 3235a8f2a6653100a81d56e630b716240673059e..a3173f07836b611fbde1b09a40665c7a21fd3fce 100644
--- a/Geo/Geo.cpp
+++ b/Geo/Geo.cpp
@@ -1178,7 +1178,6 @@ void DeleteShape(int Type, int Num)
if(gf) GModel::current()->remove(gf);
}
break;
-
case MSH_VOLUME_FROM_GMODEL:
{
GRegion *gr = GModel::current()->getRegionByTag(Num);
diff --git a/Geo/discreteEdge.cpp b/Geo/discreteEdge.cpp
index c112f4df7d47da90ab358fec1f6b76985c99794b..d94a56b2640bf623b09b5e64daa348422a0744c9 100644
--- a/Geo/discreteEdge.cpp
+++ b/Geo/discreteEdge.cpp
@@ -490,3 +490,9 @@ Range<double> discreteEdge::parBounds(int i) const
{
return Range<double>(0, lines.size());
}
+
+void discreteEdge::writeGEO(FILE *fp)
+{
+ if (getBeginVertex() && getEndVertex())
+ fprintf(fp, "Discrete Line(%d) = {%d,%d};\n",tag(), getBeginVertex()->tag(), getEndVertex()->tag());
+}
diff --git a/Geo/discreteEdge.h b/Geo/discreteEdge.h
index 5c475bc28bc051d712cdaeaf608d19e647952135..29671932c147a48a0a92b289837c01c1a7572529 100644
--- a/Geo/discreteEdge.h
+++ b/Geo/discreteEdge.h
@@ -32,6 +32,7 @@ class discreteEdge : public GEdge {
void setBoundVertices();
void createTopo();
void computeNormals () const;
+ void writeGEO(FILE *fp);
};
#endif
diff --git a/Geo/discreteFace.cpp b/Geo/discreteFace.cpp
index 5a56ea02efb1c4f85cb3ad0fa3117acc976266ae..d872593d5f8a0e1da7b6855541d34af749cd9840 100644
--- a/Geo/discreteFace.cpp
+++ b/Geo/discreteFace.cpp
@@ -122,3 +122,17 @@ void discreteFace::secondDer(const SPoint2 ¶m,
return;
}
}
+
+void discreteFace::writeGEO(FILE *fp)
+{
+ fprintf(fp, "Discrete Face(%d) = {",tag());
+ int count = 0;
+ for (std::list<GEdge*>::iterator it = l_edges.begin();
+ it != l_edges.end() ;++it){
+ if (count == 0)fprintf(fp, "%d",(*it)->tag());
+ else fprintf(fp, ",%d",(*it)->tag());
+ count ++;
+ }
+ fprintf(fp, "};\n");
+
+}
diff --git a/Geo/discreteFace.h b/Geo/discreteFace.h
index db16cfa9a41279f6f3ee0c26b5b4ad53b708d4cb..78bfa782e9ca7ae22aeb03e61f1520b471fedb31 100644
--- a/Geo/discreteFace.h
+++ b/Geo/discreteFace.h
@@ -25,6 +25,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);
+ void discreteFace::writeGEO(FILE *fp);
};
#endif
diff --git a/Mesh/Generator.cpp b/Mesh/Generator.cpp
index 3bebd69de2aed13756a788be0c3ccc8673717fa5..276f6b84703b1fea73ffb78a4e2299a8918625a7 100644
--- a/Mesh/Generator.cpp
+++ b/Mesh/Generator.cpp
@@ -616,7 +616,7 @@ void GenerateMesh(GModel *m, int ask)
}
// Subdivide into quads or hexas
- if(m->getMeshStatus() == 2 && CTX::instance()->mesh.algoSubdivide == 1)
+ if(m->getMeshStatus() == 2 && CTX::instance()->mesh.recombineAll && CTX::instance()->mesh.algoRecombine != 1)
RefineMesh(m, CTX::instance()->mesh.secondOrderLinear, true);
else if(m->getMeshStatus() == 3 && CTX::instance()->mesh.algoSubdivide == 2)
RefineMesh(m, CTX::instance()->mesh.secondOrderLinear, false, true);
diff --git a/Mesh/meshGFace.cpp b/Mesh/meshGFace.cpp
index dd8da594edd71d6a9352cce5efcc35812f71c5c0..2c6be0d2950c1c02ae33ac2f336e26b128b2b000 100644
--- a/Mesh/meshGFace.cpp
+++ b/Mesh/meshGFace.cpp
@@ -296,7 +296,7 @@ static void remeshUnrecoveredEdges(std::map<MVertex*, BDS_Point*> &recoverMapInv
static bool algoDelaunay2D(GFace *gf)
{
// HACK REMOVE ME
- if (gf->geomType() == GEntity::CompoundSurface)return true;
+ // if (gf->geomType() == GEntity::CompoundSurface)return true;
if(noSeam(gf) && (CTX::instance()->mesh.algo2d == ALGO_2D_DELAUNAY ||
CTX::instance()->mesh.algo2d == ALGO_2D_BAMG ||
@@ -872,8 +872,8 @@ static bool meshGenerator(GFace *gf, int RECUR_ITER,
// BDS mesh is passed in order not to recompute local coordinates of
// vertices
if(algoDelaunay2D(gf)){
- if(CTX::instance()->mesh.algo2d == ALGO_2D_FRONTAL ||
- gf->geomType() == GEntity::CompoundSurface)
+ if(CTX::instance()->mesh.algo2d == ALGO_2D_FRONTAL)// ||
+ // gf->geomType() == GEntity::CompoundSurface)
bowyerWatsonFrontal(gf);
else if(CTX::instance()->mesh.algo2d == ALGO_2D_DELAUNAY)
bowyerWatson(gf);
@@ -1454,8 +1454,8 @@ static bool meshGeneratorPeriodic(GFace *gf, bool debug = true)
}
if(algoDelaunay2D(gf)){
- if(CTX::instance()->mesh.algo2d == ALGO_2D_FRONTAL ||
- gf->geomType() == GEntity::CompoundSurface)
+ if(CTX::instance()->mesh.algo2d == ALGO_2D_FRONTAL)// ||
+ // gf->geomType() == GEntity::CompoundSurface)
bowyerWatsonFrontal(gf);
else if(CTX::instance()->mesh.algo2d == ALGO_2D_DELAUNAY)
bowyerWatson(gf);
diff --git a/Mesh/meshGFaceDelaunayInsertion.cpp b/Mesh/meshGFaceDelaunayInsertion.cpp
index 44452b2c8003123f334cfdb2661ad5ffadd88864..6c8bec0890cad8fd20ca448eef4533bd3d94d434 100644
--- a/Mesh/meshGFaceDelaunayInsertion.cpp
+++ b/Mesh/meshGFaceDelaunayInsertion.cpp
@@ -16,7 +16,7 @@
#include "Numeric.h"
#include "STensor3.h"
-const double LIMIT_ = 0.5 * sqrt(2.0);
+const double LIMIT_ = 0.5 * sqrt(2.0) * 1;
static const bool _experimental_anisotropic_blues_band_ = false;
diff --git a/Mesh/meshGFaceOptimize.cpp b/Mesh/meshGFaceOptimize.cpp
index dbf4610d1e322fa9de5a0cacb0c5a38da5c2fc4c..cea0fc23a2812c9cc9b13a260f65127a38b610f8 100644
--- a/Mesh/meshGFaceOptimize.cpp
+++ b/Mesh/meshGFaceOptimize.cpp
@@ -849,6 +849,10 @@ void laplaceSmoothing(GFace *gf)
int _edgeSwapQuadsForBetterQuality ( GFace *gf ) {
return 0;
+
+ v2t_cont adj_v;
+ buildVertexToElement(gf->quadrangles,adj_v);
+
e2t_cont adj;
// buildEdgeToElement(gf->triangles, adj);
buildEdgeToElement(gf->quadrangles, adj);
@@ -882,30 +886,41 @@ int _edgeSwapQuadsForBetterQuality ( GFace *gf ) {
if (ed.getVertex(0) == v2 && ed.getVertex(1) != v1)v22 = ed.getVertex(1);
if (ed.getVertex(1) == v2 && ed.getVertex(0) != v1)v22 = ed.getVertex(0);
}
- MQuadrangle *q1 = new MQuadrangle (v11,v22,v2,v12);
- MQuadrangle *q2 = new MQuadrangle (v22,v11,v1,v21);
-
- double sold = surfaceFaceUV(e1,gf) + surfaceFaceUV(e2,gf);
- double snew = surfaceFaceUV(q1,gf) + surfaceFaceUV(q2,gf);
-
- double min_quality_old = std::min(e1->etaShapeMeasure(),e2->etaShapeMeasure());
- double min_quality_new = std::min(q1->etaShapeMeasure(),q2->etaShapeMeasure());
-
- // if (min_quality_old < min_quality_new){
- if (sold > 1.00001*snew){
- printf("%g %g\n",sold, snew);
+ MQuadrangle *q1=0,*q2=0;
+ v2t_cont :: iterator it2 = adj_v.find(v2);
+ v2t_cont :: iterator it1 = adj_v.find(v1);
+ v2t_cont :: iterator it22 = adj_v.find(v22);
+ v2t_cont :: iterator it12 = adj_v.find(v12);
+ v2t_cont :: iterator it21 = adj_v.find(v21);
+ v2t_cont :: iterator it11 = adj_v.find(v11);
+ /*
+ int connectivity_default_config_1 =
+ abs(it1->size() - 4) +
+ abs(it2->size() - 4) +
+ abs(it2->size() - 4) +
+ abs(it2->size() - 4) ;
+ */
+ if (it2->second.size() >= 5 && it1->second.size() >= 5){
+ if (it22->second.size() <= 3 && it11->second.size() <= 3) {
+ q1 = new MQuadrangle (v11,v22,v2,v12);
+ q2 = new MQuadrangle (v22,v11,v1,v21);
+ }
+ else if (it21->second.size() <= 3 && it12->second.size() <= 3) {
+ q1 = new MQuadrangle (v11,v12,v21,v1);
+ q2 = new MQuadrangle (v21,v12,v2,v22);
+ }
+ }
+ if (q1 && q2){
deleted.insert(e1);
deleted.insert(e2);
created.push_back(q1);
created.push_back(q2);
+ printf("edge swap performed\n");
COUNT++;
}
- else{
- delete q1;
- delete q2;
- }
}
}
+ if (COUNT)break;
}
for(unsigned int i = 0; i < gf->quadrangles.size(); i++){
@@ -1371,6 +1386,19 @@ static int _recombineIntoQuads(GFace *gf, int recur_level, bool cubicGraph = 1)
elist[2*i] = t2n[pairs[i].t1];
elist[2*i+1] = t2n[pairs[i].t2];
elen [i] = (int) pairs[i].angle;
+ double angle = atan2(pairs[i].n1->y()-pairs[i].n2->y(),
+ pairs[i].n1->x()-pairs[i].n2->x());
+
+ double x = .5*(pairs[i].n1->x()+pairs[i].n2->x());
+ double y = .5*(pairs[i].n1->y()+pairs[i].n2->y());
+ double opti = atan2(y,x);
+ angle -= opti;
+ while (angle < 0 || angle > M_PI/2){
+ if (angle < 0) angle += M_PI/2;
+ if (angle > M_PI/2) angle -= M_PI/2;
+ }
+ //elen [i] = (int) 180. * fabs(angle-M_PI/4)/M_PI;
+
int NB = 0;
if (pairs[i].n1->onWhat()->dim() < 2)NB++;
if (pairs[i].n2->onWhat()->dim() < 2)NB++;
@@ -1558,7 +1586,7 @@ static int _recombineIntoQuads(GFace *gf, int recur_level, bool cubicGraph = 1)
void recombineIntoQuads(GFace *gf)
{
- // gf->model()->writeMSH("before.msh");
+ gf->model()->writeMSH("before.msh");
removeFourTrianglesNodes(gf, false);
int success = _recombineIntoQuads(gf,0);
gf->model()->writeMSH("raw.msh");
diff --git a/Mesh/meshGFaceOptimize.h b/Mesh/meshGFaceOptimize.h
index e6b8e7538a5ee6a848da596ad29d24abc3688f92..bda526530abb7e3841b51e3bafd82448aa60a1e6 100644
--- a/Mesh/meshGFaceOptimize.h
+++ b/Mesh/meshGFaceOptimize.h
@@ -30,6 +30,9 @@ class edge_angle {
typedef std::map<MVertex*, std::vector<MElement*> > v2t_cont;
typedef std::map<MEdge, std::pair<MElement*, MElement*>, Less_Edge> e2t_cont;
+template <class T> void buildVertexToElement(std::vector<T*> &eles, v2t_cont &adj);
+template <class T> void buildEdgeToElement(std::vector<T*> &eles, e2t_cont &adj);
+
void buildVertexToTriangle(std::vector<MTriangle*> &, v2t_cont &adj);
void buildEdgeToTriangle(std::vector<MTriangle*> &, e2t_cont &adj);
void buildListOfEdgeAngle(e2t_cont adj, std::vector<edge_angle> &edges_detected,
diff --git a/Mesh/meshRefine.cpp b/Mesh/meshRefine.cpp
index c2e526392b43ca26891757af1834e5728372b444..021c17ca9a7ff6be31a5f32478f0fe53723cf07c 100644
--- a/Mesh/meshRefine.cpp
+++ b/Mesh/meshRefine.cpp
@@ -17,6 +17,8 @@
#include "MPyramid.h"
#include "GmshMessage.h"
#include "OS.h"
+#include "Context.h"
+#include "meshGFaceOptimize.h"
class MVertexLessThanParam{
public:
@@ -375,8 +377,11 @@ void RefineMesh(GModel *m, bool linear, bool splitIntoQuads, bool splitIntoHexas
// mesh
for(GModel::eiter it = m->firstEdge(); it != m->lastEdge(); ++it)
Subdivide(*it);
- for(GModel::fiter it = m->firstFace(); it != m->lastFace(); ++it)
+ for(GModel::fiter it = m->firstFace(); it != m->lastFace(); ++it){
Subdivide(*it, splitIntoQuads, splitIntoHexas, faceVertices);
+ for(int i = 0; i < CTX::instance()->mesh.nbSmoothing; i++)
+ laplaceSmoothing(*it);
+ }
for(GModel::riter it = m->firstRegion(); it != m->lastRegion(); ++it)
Subdivide(*it, splitIntoHexas, faceVertices);
diff --git a/benchmarks/boolean/wikipedia.lua b/benchmarks/boolean/wikipedia.lua
index bd1ee4e9b56f403e2ca28a3d1162b6e613dd7d15..1d7a612cd5cc8fca0329e8a61ba5efcac19dc9b7 100644
--- a/benchmarks/boolean/wikipedia.lua
+++ b/benchmarks/boolean/wikipedia.lua
@@ -1,8 +1,8 @@
-- from http://en.wikipedia.org/wiki/Constructive_solid_geometry
-R = 1.0;
-s = .8;
+R = 1.4;
+s = .4;
t = 1.35;
myModel = GModel();
myModel:addBlock({-R,-R,-R},{R,R,R});