From 6ab4d46a8437be1e32eb675cde587b5e482fd8a7 Mon Sep 17 00:00:00 2001
From: Jean-Francois Remacle <jean-francois.remacle@uclouvain.be>
Date: Wed, 3 Nov 2010 16:24:20 +0000
Subject: [PATCH] added hessian based error estimation & adaptation better quad
meshing
---
Common/OctreeInternals.cpp | 2 +-
Fltk/statisticsWindow.cpp | 2 -
Geo/GModel.cpp | 1 -
Geo/MElement.cpp | 12 +
Geo/MElement.h | 1 +
Geo/MElementOctree.cpp | 28 ++-
Geo/MElementOctree.h | 1 +
Geo/MQuadrangle.cpp | 12 +
Geo/OCCEdge.cpp | 2 +-
Mesh/BDS.cpp | 6 +-
Mesh/BDS.h | 2 +-
Mesh/BackgroundMesh.cpp | 5 +-
Mesh/meshGFace.cpp | 28 ++-
Mesh/meshGFaceOptimize.cpp | 391 +++++++++++++++++++----------
Mesh/meshGFaceOptimize.h | 23 +-
Mesh/meshGFaceQuadrilateralize.cpp | 3 +-
Mesh/meshRefine.cpp | 6 +-
benchmarks/boolean/sphere.lua | 4 +-
18 files changed, 371 insertions(+), 158 deletions(-)
diff --git a/Common/OctreeInternals.cpp b/Common/OctreeInternals.cpp
index a12113f9a2..2c500d379a 100644
--- a/Common/OctreeInternals.cpp
+++ b/Common/OctreeInternals.cpp
@@ -417,6 +417,6 @@ void *searchAllElements(octantBucket *_buckets_head, double *_pt, globalInfo *_g
if (flag1)
return (void *)(_elements);
- Msg::Warning("This point is not found in any element! It is not in the domain");
+ // Msg::Warning("This point is not found in any element! It is not in the domain");
return NULL;
}
diff --git a/Fltk/statisticsWindow.cpp b/Fltk/statisticsWindow.cpp
index beeb11201c..5a431cc4e0 100644
--- a/Fltk/statisticsWindow.cpp
+++ b/Fltk/statisticsWindow.cpp
@@ -205,8 +205,6 @@ 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;
diff --git a/Geo/GModel.cpp b/Geo/GModel.cpp
index d63236b4ca..9697d63b11 100644
--- a/Geo/GModel.cpp
+++ b/Geo/GModel.cpp
@@ -2508,5 +2508,4 @@ void GModel::registerBindings(binding *b)
cm->setDescription("Assigns partition tags to boundary elements. Should be called "
"only after the partitions have been assigned");
cm->setArgNames("createGhostCells",NULL);
-
}
diff --git a/Geo/MElement.cpp b/Geo/MElement.cpp
index d8416bafe5..1d36d466a9 100644
--- a/Geo/MElement.cpp
+++ b/Geo/MElement.cpp
@@ -1095,6 +1095,15 @@ const fullMatrix<double> &MElement::getGradShapeFunctionsAtNodes (int functionSp
return mat;
}
+void MElement::xyzTouvw(fullMatrix<double> *xu){
+ double _xyz[3] = {(*xu)(0,0),(*xu)(0,1),(*xu)(0,2)},_uvw[3] ;
+ xyz2uvw(_xyz,_uvw);
+ (*xu)(1,0) = _uvw[0];
+ (*xu)(1,1) = _uvw[1];
+ (*xu)(1,2) = _uvw[2];
+}
+
+
#include "Bindings.h"
void MElement::registerBindings(binding *b)
@@ -1125,4 +1134,7 @@ void MElement::registerBindings(binding *b)
cm = cb->addMethod("getJacobianDeterminant", &MElement::getJacobianDeterminant);
cm->setDescription("return the jacobian of the determinant of the transformation");
cm->setArgNames("u","v","w",NULL);
+ cm = cb->addMethod("xyzTouvw", &MElement::xyzTouvw);
+ cm->setDescription("get uvw from xyz");
+ cm->setArgNames("xyzuvw",NULL);
}
diff --git a/Geo/MElement.h b/Geo/MElement.h
index f758c8c166..5f7a974ac2 100644
--- a/Geo/MElement.h
+++ b/Geo/MElement.h
@@ -244,6 +244,7 @@ class MElement
// invert the parametrisation
virtual void xyz2uvw(double xyz[3], double uvw[3]);
+ void xyzTouvw(fullMatrix<double> *xu);
// move point between parent and element parametric spaces
virtual void movePointFromParentSpaceToElementSpace(double &u, double &v, double &w);
diff --git a/Geo/MElementOctree.cpp b/Geo/MElementOctree.cpp
index c74568b713..9e2d94aa50 100644
--- a/Geo/MElementOctree.cpp
+++ b/Geo/MElementOctree.cpp
@@ -55,7 +55,7 @@ static int MElementInEle(void *a, double *x)
return e->isInside(uvw[0], uvw[1], uvw[2]) ? 1 : 0;
}
-MElementOctree::MElementOctree(GModel *m)
+MElementOctree::MElementOctree(GModel *m) : _gm(m)
{
SBoundingBox3d bb = m->bounds();
// make bounding box larger up to (absolute) geometrical tolerance
@@ -80,9 +80,11 @@ MElementOctree::MElementOctree(GModel *m)
MElementOctree::MElementOctree(std::vector<MElement*> &v)
{
+ _gm = 0;
SBoundingBox3d bb;
for (unsigned int i = 0; i < v.size(); i++){
for(int j = 0; j < v[i]->getNumVertices(); j++){
+ if (!_gm)_gm = v[i]->getVertex(j)->onWhat()->model();
bb += SPoint3(v[i]->getVertex(j)->x(),
v[i]->getVertex(j)->y(),
v[i]->getVertex(j)->z());
@@ -114,6 +116,30 @@ MElement *MElementOctree::find(double x, double y, double z, int dim)
{
double P[3] = {x, y, z};
MElement *e = (MElement*)Octree_Search(P, _octree);
+
+ if (!e){
+ double initialTol = 1.e-6;
+ double tol = initialTol;
+ while ( tol < .1){
+ tol *= 10;
+ MElement::setTolerance(tol);
+ std::vector<GEntity*> entities;
+ _gm->getEntities(entities);
+ for(unsigned int i = 0; i < entities.size(); i++){
+ for(unsigned int j = 0; j < entities[i]->getNumMeshElements(); j++){
+ e = entities[i]->getMeshElement(j);
+ if (dim == -1 || e->getDim() == dim){
+ if (MElementInEle(e, P)){
+ MElement::setTolerance(initialTol);
+ return e;
+ }
+ }
+ }
+ }
+ }
+ }
+
+
if (dim == -1 || !e || e->getDim() == dim)
return e;
std::list<void*> l;
diff --git a/Geo/MElementOctree.h b/Geo/MElementOctree.h
index 112589611f..4b8511f52f 100644
--- a/Geo/MElementOctree.h
+++ b/Geo/MElementOctree.h
@@ -15,6 +15,7 @@ class MElement;
class MElementOctree{
private:
Octree *_octree;
+ GModel *_gm;
public:
MElementOctree(GModel *);
MElementOctree(std::vector<MElement*> &);
diff --git a/Geo/MQuadrangle.cpp b/Geo/MQuadrangle.cpp
index 0563ae6d29..2d74ebfc12 100644
--- a/Geo/MQuadrangle.cpp
+++ b/Geo/MQuadrangle.cpp
@@ -219,6 +219,18 @@ void MQuadrangle::getIntegrationPoints(int pOrder, int *npts, IntPt **pts)
double MQuadrangle::etaShapeMeasure()
{
+ SVector3 v01 (_v[1]->x()-_v[0]->x(),_v[1]->y()-_v[0]->y(),_v[1]->z()-_v[0]->z());
+ SVector3 v12 (_v[2]->x()-_v[1]->x(),_v[2]->y()-_v[1]->y(),_v[2]->z()-_v[1]->z());
+ SVector3 v23 (_v[3]->x()-_v[2]->x(),_v[3]->y()-_v[2]->y(),_v[3]->z()-_v[2]->z());
+ SVector3 v30 (_v[0]->x()-_v[3]->x(),_v[0]->y()-_v[3]->y(),_v[0]->z()-_v[3]->z());
+
+ SVector3 a = crossprod(v01,v12);
+ SVector3 b = crossprod(v12,v23);
+ SVector3 c = crossprod(v23,v30);
+ SVector3 d = crossprod(v30,v01);
+
+ if (dot(a,b) < 0 || dot(a,c) < 0 || dot(a,d) < 0 )return 0.0;
+
double a1 = 180 * angle3Vertices(_v[0], _v[1], _v[2]) / M_PI;
double a2 = 180 * angle3Vertices(_v[1], _v[2], _v[3]) / M_PI;
double a3 = 180 * angle3Vertices(_v[2], _v[3], _v[0]) / M_PI;
diff --git a/Geo/OCCEdge.cpp b/Geo/OCCEdge.cpp
index 15b14e8196..142413f49f 100644
--- a/Geo/OCCEdge.cpp
+++ b/Geo/OCCEdge.cpp
@@ -97,7 +97,7 @@ SPoint2 OCCEdge::reparamOnFace(const GFace *face, double epar, int dir) const
const double dx = p1.x()-p2.x();
const double dy = p1.y()-p2.y();
const double dz = p1.z()-p2.z();
- if(sqrt(dx * dx + dy * dy + dz * dz) > 1.e-4 * CTX::instance()->lc){
+ if(sqrt(dx * dx + dy * dy + dz * dz) > 1.e-2 * CTX::instance()->lc && 0){
// return reparamOnFace(face, epar,-1);
Msg::Warning("Reparam on face partially failed for curve %d surface %d at point %g",
tag(), face->tag(), epar);
diff --git a/Mesh/BDS.cpp b/Mesh/BDS.cpp
index 1c98232a95..d25cb34423 100644
--- a/Mesh/BDS.cpp
+++ b/Mesh/BDS.cpp
@@ -263,10 +263,11 @@ BDS_Edge *BDS_Mesh::recover_edge_fast(BDS_Point *p1, BDS_Point *p2){
return 0;
}
-BDS_Edge *BDS_Mesh::recover_edge(int num1, int num2, std::set<EdgeToRecover> *e2r,
+BDS_Edge *BDS_Mesh::recover_edge(int num1, int num2, bool &_fatal, std::set<EdgeToRecover> *e2r,
std::set<EdgeToRecover> *not_recovered)
{
BDS_Edge *e = find_edge(num1, num2);
+ _fatal = false;
if(e) return e;
@@ -331,6 +332,7 @@ BDS_Edge *BDS_Mesh::recover_edge(int num1, int num2, std::set<EdgeToRecover> *e2
outputScalarField(triangles, "debugr.pos", 0);
Msg::Debug("edge %d %d cannot be recovered at all, look at debugp.pos "
"and debugr.pos", num1, num2);
+ _fatal = true;
return 0;
}
return eee;
@@ -1170,7 +1172,7 @@ bool BDS_Mesh::smooth_point_centroid(BDS_Point *p, GFace *gf, bool test_quality)
}
/* TEST */
- double R; SPoint3 c; bool isSphere = gf->isSphere(R,c);
+ bool isSphere = gf->geomType() == GEntity::Sphere;
double XX=0,YY=0,ZZ=0;
double U = 0;
diff --git a/Mesh/BDS.h b/Mesh/BDS.h
index 937a75d9ef..4542f69c1f 100644
--- a/Mesh/BDS.h
+++ b/Mesh/BDS.h
@@ -440,7 +440,7 @@ class BDS_Mesh
void add_geom(int degree, int tag);
BDS_GeomEntity *get_geom(int p1, int p2);
// 2D operators
- BDS_Edge *recover_edge(int p1, int p2, std::set<EdgeToRecover> *e2r=0,
+ BDS_Edge *recover_edge(int p1, int p2, bool &_fatal, std::set<EdgeToRecover> *e2r=0,
std::set<EdgeToRecover> *not_recovered=0);
BDS_Edge *recover_edge_fast(BDS_Point *p1, BDS_Point *p2);
bool swap_edge(BDS_Edge*, const BDS_SwapEdgeTest &theTest);
diff --git a/Mesh/BackgroundMesh.cpp b/Mesh/BackgroundMesh.cpp
index aeaf804e12..026f80a21f 100644
--- a/Mesh/BackgroundMesh.cpp
+++ b/Mesh/BackgroundMesh.cpp
@@ -250,9 +250,10 @@ double BGM_MeshSize(GEntity *ge, double U, double V,
// lc from fields
double l4 = MAX_LC;
- FieldManager *fields = GModel::current()->getFields();
+ FieldManager *fields = ge->model()->getFields();
if(fields->background_field > 0){
Field *f = fields->get(fields->background_field);
+ // printf("field %p %s %d %p\n",f,f->getName(),fields->size(), ge->model());
if(f) l4 = (*f)(X, Y, Z, ge);
}
@@ -294,7 +295,7 @@ SMetric3 BGM_MeshMetric(GEntity *ge,
// lc from fields
SMetric3 l4(1./(MAX_LC*MAX_LC));
- FieldManager *fields = GModel::current()->getFields();
+ FieldManager *fields = ge->model()->getFields();
if(fields->background_field > 0){
Field *f = fields->get(fields->background_field);
if(f){
diff --git a/Mesh/meshGFace.cpp b/Mesh/meshGFace.cpp
index a1737ba6fc..a4b0c830af 100644
--- a/Mesh/meshGFace.cpp
+++ b/Mesh/meshGFace.cpp
@@ -335,6 +335,8 @@ static bool recoverEdge(BDS_Mesh *m, GEdge *ge,
g = m->get_geom(ge->tag(), 1);
}
+ bool _fatallyFailed;
+
for(unsigned int i = 0; i < ge->lines.size(); i++){
MVertex *vstart = ge->lines[i]->getVertex(0);
MVertex *vend = ge->lines[i]->getVertex(1);
@@ -346,14 +348,14 @@ static bool recoverEdge(BDS_Mesh *m, GEdge *ge,
if(pass == 1)
e2r->insert(EdgeToRecover(pstart->iD, pend->iD, ge));
else{
- BDS_Edge *e = m->recover_edge(pstart->iD, pend->iD, e2r, notRecovered);
+ BDS_Edge *e = m->recover_edge(pstart->iD, pend->iD, _fatallyFailed, e2r, notRecovered);
if(e) e->g = g;
- //else {
- // Msg::Error("Unable to recover an edge %g %g && %g %g (%d/%d)",
- // vstart->x(), vstart->y(), vend->x(), vend->y(), i,
- // ge->mesh_vertices.size());
- // return false;
- // }
+ else {
+ if (_fatallyFailed) Msg::Error("Unable to recover an edge %g %g && %g %g (%d/%d)",
+ vstart->x(), vstart->y(), vend->x(), vend->y(), i,
+ ge->mesh_vertices.size());
+ return !_fatallyFailed;
+ }
}
}
}
@@ -575,7 +577,11 @@ static bool meshGenerator(GFace *gf, int RECUR_ITER,
ite = edges.begin();
while(ite != edges.end()){
if(!(*ite)->isMeshDegenerated()){
- recoverEdge(m, *ite, recoverMapInv, &edgesToRecover, &edgesNotRecovered, 2);
+ if (!recoverEdge(m, *ite, recoverMapInv, &edgesToRecover, &edgesNotRecovered, 2)){
+ delete m;
+ gf->meshStatistics.status = GFace::FAILED;
+ return false;
+ }
}
++ite;
}
@@ -900,7 +906,7 @@ static bool meshGenerator(GFace *gf, int RECUR_ITER,
if((CTX::instance()->mesh.recombineAll || gf->meshAttributes.recombine) &&
!CTX::instance()->mesh.optimizeLloyd)
- recombineIntoQuads(gf);
+ recombineIntoQuadsIterative(gf);
computeElementShapes(gf, gf->meshStatistics.worst_element_shape,
gf->meshStatistics.average_element_shape,
@@ -1283,11 +1289,13 @@ static bool meshGeneratorPeriodic(GFace *gf, bool debug = true)
outputScalarField(m->triangles, name, 1);
}
+ bool _fatallyFailed;
+
for(unsigned int i = 0; i < edgeLoops_BDS.size(); i++){
std::vector<BDS_Point*> &edgeLoop_BDS = edgeLoops_BDS[i];
for(unsigned int j = 0; j < edgeLoop_BDS.size(); j++){
BDS_Edge * e = m->recover_edge
- (edgeLoop_BDS[j]->iD, edgeLoop_BDS[(j + 1) % edgeLoop_BDS.size()]->iD);
+ (edgeLoop_BDS[j]->iD, edgeLoop_BDS[(j + 1) % edgeLoop_BDS.size()]->iD, _fatallyFailed);
if(!e){
Msg::Error("Impossible to recover the edge %d %d", edgeLoop_BDS[j]->iD,
edgeLoop_BDS[(j + 1) % edgeLoop_BDS.size()]->iD);
diff --git a/Mesh/meshGFaceOptimize.cpp b/Mesh/meshGFaceOptimize.cpp
index b6f5180171..f3be29529e 100644
--- a/Mesh/meshGFaceOptimize.cpp
+++ b/Mesh/meshGFaceOptimize.cpp
@@ -187,26 +187,6 @@ void transferDataStructure(GFace *gf, std::set<MTri3*, compareTri3Ptr> &AllTris,
}
}
-template <class T>
-void buildVertexToElement(std::vector<T*> &eles, v2t_cont &adj)
-{
- for(unsigned int i = 0; i < eles.size(); i++){
- T *t = eles[i];
- for(int j = 0; j < t->getNumVertices(); j++){
- MVertex *v = t->getVertex(j);
- v2t_cont :: iterator it = adj.find(v);
- if(it == adj.end()){
- std::vector<MElement*> one;
- one.push_back(t);
- adj[v] = one;
- }
- else{
- it->second.push_back(t);
- }
- }
- }
-}
-
void buildVertexToTriangle(std::vector<MTriangle*> &eles, v2t_cont &adj)
{
adj.clear();
@@ -569,7 +549,7 @@ static bool _isItAGoodIdeaToCollapseThatVertex (GFace *gf,
static bool _isItAGoodIdeaToMoveThatVertex (GFace *gf,
- std::vector<MElement*> &e1,
+ const std::vector<MElement*> &e1,
MVertex *v1,
const SPoint2 &before,
const SPoint2 &after){
@@ -775,67 +755,65 @@ int removeDiamonds(GFace *gf){
return nbRemove;
}
-
+static void _relocateVertex (GFace *gf,
+ MVertex *ver,
+ const std::vector<MElement*> <){
+ double R; SPoint3 c; bool isSphere = gf->isSphere(R,c);
+ if( ver->onWhat()->dim() == 2){
+ double initu,initv;
+ ver->getParameter(0, initu);
+ ver->getParameter(1, initv);
+ double cu = 0, cv = 0;
+ double XX=0,YY=0,ZZ=0;
+ double pu[4], pv[4];
+ double fact = 0.0;
+ for(unsigned int i = 0; i < lt.size(); i++){
+ parametricCoordinates(lt[i], gf, pu, pv, ver);
+ cu += (pu[0] + pu[1] + pu[2]);
+ cv += (pv[0] + pv[1] + pv[2]);
+ XX += lt[i]->getVertex(0)->x()+lt[i]->getVertex(1)->x()+lt[i]->getVertex(2)->x();
+ YY += lt[i]->getVertex(0)->y()+lt[i]->getVertex(1)->y()+lt[i]->getVertex(2)->y();
+ ZZ += lt[i]->getVertex(0)->z()+lt[i]->getVertex(1)->z()+lt[i]->getVertex(2)->z();
+ if(lt[i]->getNumVertices() == 4){
+ cu += pu[3];
+ cv += pv[3];
+ XX += lt[i]->getVertex(3)->x();
+ YY += lt[i]->getVertex(3)->y();
+ ZZ += lt[i]->getVertex(3)->z();
+ }
+ fact += lt[i]->getNumVertices();
+ }
+ if(fact != 0.0){
+ SPoint2 before(initu,initv);
+ SPoint2 after(cu / fact,cv / fact);
+ if (isSphere){
+ GPoint gp = gf->closestPoint(SPoint3(XX/fact, YY/fact, ZZ/fact), after);
+ after = SPoint2(gp.u(),gp.v());
+ }
+ bool success = _isItAGoodIdeaToMoveThatVertex (gf, lt, ver,before,after);
+ if (success){
+ ver->setParameter(0, after.x());
+ ver->setParameter(1, after.y());
+ GPoint pt = gf->point(after);
+ if(pt.succeeded()){
+ ver->x() = pt.x();
+ ver->y() = pt.y();
+ ver->z() = pt.z();
+ }
+ }
+ }
+ }
+}
+
void laplaceSmoothing(GFace *gf)
{
v2t_cont adj;
buildVertexToElement(gf->triangles, adj);
buildVertexToElement(gf->quadrangles, adj);
-
- /* TEST */
- double R; SPoint3 c; bool isSphere = gf->isSphere(R,c);
-
for(int i = 0; i < 5; i++){
v2t_cont :: iterator it = adj.begin();
while (it != adj.end()){
- MVertex *ver= it->first;
- GEntity *ge = ver->onWhat();
- // this vertex in internal to the face
- if(ge->dim() == 2){
- double initu,initv;
- ver->getParameter(0, initu);
- ver->getParameter(1, initv);
- const std::vector<MElement*> < = it->second;
- double cu = 0, cv = 0;
- double XX=0,YY=0,ZZ=0;
- double pu[4], pv[4];
- double fact = 0.0;
- for(unsigned int i = 0; i < lt.size(); i++){
- parametricCoordinates(lt[i], gf, pu, pv, ver);
- cu += (pu[0] + pu[1] + pu[2]);
- cv += (pv[0] + pv[1] + pv[2]);
- XX += lt[i]->getVertex(0)->x()+lt[i]->getVertex(1)->x()+lt[i]->getVertex(2)->x();
- YY += lt[i]->getVertex(0)->y()+lt[i]->getVertex(1)->y()+lt[i]->getVertex(2)->y();
- ZZ += lt[i]->getVertex(0)->z()+lt[i]->getVertex(1)->z()+lt[i]->getVertex(2)->z();
- if(lt[i]->getNumVertices() == 4){
- cu += pu[3];
- cv += pv[3];
- XX += lt[i]->getVertex(3)->x();
- YY += lt[i]->getVertex(3)->y();
- ZZ += lt[i]->getVertex(3)->z();
- }
- fact += lt[i]->getNumVertices();
- }
- if(fact != 0.0){
- SPoint2 before(initu,initv);
- SPoint2 after(cu / fact,cv / fact);
- if (isSphere){
- GPoint gp = gf->closestPoint(SPoint3(XX/fact, YY/fact, ZZ/fact), after);
- after = SPoint2(gp.u(),gp.v());
- }
- bool success = _isItAGoodIdeaToMoveThatVertex (gf, it->second, ver,before,after);
- if (success){
- ver->setParameter(0, after.x());
- ver->setParameter(1, after.y());
- GPoint pt = gf->point(after);
- if(pt.succeeded()){
- ver->x() = pt.x();
- ver->y() = pt.y();
- ver->z() = pt.z();
- }
- }
- }
- }
+ _relocateVertex(gf,it->first,it->second);
++it;
}
}
@@ -843,16 +821,11 @@ 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);
- std::set<MElement*>touched;
std::vector<MQuadrangle*>created;
std::set<MElement*>deleted;
@@ -866,8 +839,12 @@ int _edgeSwapQuadsForBetterQuality ( GFace *gf ) {
MVertex *v2 = it->first.getVertex(1);
MElement *e1 = it->second.first;
MElement *e2 = it->second.second;
- if (deleted.find(e1) == deleted.end() ||
- deleted.find(e2) == deleted.end()){
+
+ double worst_quality_old = std::min(e1->etaShapeMeasure(),e2->etaShapeMeasure());
+
+ if (worst_quality_old < .1 && (
+ deleted.find(e1) == deleted.end() ||
+ deleted.find(e2) == deleted.end())){
MVertex *v12,*v11,*v22,*v21;
for (int i=0;i<4;i++){
MEdge ed = e1->getEdge(i);
@@ -881,38 +858,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=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);
- }
+
+ MQuadrangle *q1A = new MQuadrangle (v11,v22,v2,v12);
+ MQuadrangle *q2A = new MQuadrangle (v22,v11,v1,v21);
+ MQuadrangle *q1B = new MQuadrangle (v11,v12,v21,v1);
+ MQuadrangle *q2B = new MQuadrangle (v21,v12,v2,v22);
+ double worst_quality_A = std::min(q1A->etaShapeMeasure(),q2A->etaShapeMeasure());
+ double worst_quality_B = std::min(q1B->etaShapeMeasure(),q2B->etaShapeMeasure());
+ // printf("worst_quality_old = %g worst_quality_A = %g worst_quality_B = %g\n",worst_quality_old,worst_quality_A,worst_quality_B);
+
+ if (0.8*worst_quality_A > worst_quality_old && 0.8*worst_quality_A > worst_quality_B){
+ deleted.insert(e1);
+ deleted.insert(e2);
+ created.push_back(q1A);
+ created.push_back(q2A);
+ delete q1B;
+ delete q2B;
+ // printf("edge swap performed\n");
+ COUNT++;
}
- if (q1 && q2){
+ else if (0.8*worst_quality_B > worst_quality_old && 0.8*worst_quality_B > worst_quality_A){
deleted.insert(e1);
deleted.insert(e2);
- created.push_back(q1);
- created.push_back(q2);
- printf("edge swap performed\n");
+ created.push_back(q1B);
+ created.push_back(q2B);
+ delete q1A;
+ delete q2A;
+ // printf("edge swap performed\n");
COUNT++;
}
+ else {
+ delete q1A;
+ delete q2A;
+ delete q1B;
+ delete q2B;
+ }
}
}
if (COUNT)break;
@@ -927,10 +907,116 @@ int _edgeSwapQuadsForBetterQuality ( GFace *gf ) {
}
}
gf->quadrangles = created;
- Msg::Debug("%i swap quads performed",COUNT);
return COUNT;
}
+static int edgeSwapQuadsForBetterQuality ( GFace *gf ) {
+ return 0;
+ int COUNT = 0;
+ while(1){
+ int k = _edgeSwapQuadsForBetterQuality (gf);
+ COUNT += k;
+ if (!k || COUNT > 10)break;
+ }
+ Msg::Debug("%i swap quads performed",COUNT);
+}
+
+
+
+int postProcessExtraEdges (GFace *gf, std::vector<std::pair<MElement*,MElement*> > &toProcess){
+
+ /*
+ v2t_cont adj;
+ buildVertexToElement(gf->triangles, adj);
+ buildVertexToElement(gf->quadrangles, adj);
+
+ std::set<MElement*>deleted;
+
+ for (int i=0; i<toProcess.size() ; i++){
+ MElement *t1 = toProcess[i].first;
+ MElement *t2 = toProcess[i].second;
+ MVertex *common = 0;
+ for (int j=0;j<3;j++){
+ if (t1->getVertex(j) == t2->getVertex(0) ||
+ t1->getVertex(j) == t2->getVertex(1) ||
+ t1->getVertex(j) == t2->getVertex(2) ){
+ common = t1->getVertex(j);
+ break;
+ }
+ }
+ if (common){
+
+ deleted.insert(t1);
+ deleted.insert(t2);
+ v2t_cont :: iterator it = adj.find(common);
+ if(it != adj.end()){
+ SPoint2 p;
+ bool success = reparamMeshVertexOnFace(it->first,gf, p);
+ MFaceVertex *newVertex = new MVertex (it->first->x(),
+ it->first->y(),
+ it->first->z(),
+ gf,
+ p.x(),
+ p.y());
+ int start1,start2;
+ int orientation1, orientation2;
+ for (int k=0;k<3;k++){
+ if (t1->getVertex(k) == it->first)start1 = k;
+ if (t2->getVertex(k) == it->first)start2 = k;
+ }
+ MQuadrangle *q1;
+ if (orientation1)
+ q1 = new MQuadrangle(newVertex,
+ t1->getVertex(start1),
+ t1->getVertex((start1+1)%3),
+ t1->getVertex((start1+2)%3));
+ else
+ q1 = new MQuadrangle(newVertex,
+ t1->getVertex((start1+2)%3),
+ t1->getVertex((start1+1)%3),
+ t1->getVertex(start1));
+
+ MQuadrangle *q2;
+ if (orientation2)
+ q2 = new MQuadrangle(newVertex,
+ t2->getVertex(start2),
+ t2->getVertex((start2+1)%3),
+ t2->getVertex((start2+2)%3));
+ else
+ q2 = new MQuadrangle(newVertex,
+ t2->getVertex((start2+2)%3),
+ t2->getVertex((start2+1)%3),
+ t2->getVertex(start2));
+
+
+ MQuadrangle *q2 = new MQuadrangle();
+ std::vector<MElement*> newAdj;
+ newAdj.push_back(q1);
+ newAdj.push_back(q2);
+ for (int k=0;k<it->second.size();k++){
+ if (it->second[k]->getNumVertices() == 4)
+ newAdj.push_back(it->second[k]);
+ }
+ gf->quadrangles.insert(q1);
+ gf->quadrangles.insert(q2);
+ _relocateVertex(gf,newVertex,newAdj);
+ }
+ }
+ }
+
+ std::vector<MTriangle*>remained;
+ for(unsigned int i = 0; i < gf->triangles.size(); i++){
+ if(deleted.find(gf->triangles[i]) == deleted.end()){
+ remained.push_back(gf->triangles[i]);
+ }
+ else {
+ delete gf->triangles[i];
+ }
+ }
+ gf->triangles = remained;
+ */
+}
+
bool edgeSwap(std::set<swapquad> &configs, MTri3 *t1, GFace *gf, int iLocalEdge,
std::vector<MTri3*> &newTris, const swapCriterion &cr,
const std::vector<double> &Us, const std::vector<double> &Vs,
@@ -1387,7 +1473,7 @@ static int _recombineIntoQuads(GFace *gf, int recur_level, bool cubicGraph = 1)
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;
+ //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;
@@ -1483,9 +1569,12 @@ static int _recombineIntoQuads(GFace *gf, int recur_level, bool cubicGraph = 1)
// FILE *f = fopen (MATCHFILE,"r");
// int nn,ne;
// fscanf(f,"%d %d",&nn,&ne);
+ std::vector<std::pair<MElement*,MElement*> > toProcess;
for (int k=0;k<elist[0];k++){
int i1 = elist[1+3*k],i2 = elist[1+3*k+1],an=elist[1+3*k+2];
- if (an == 100000){
+ // FIXME !!
+ if (an == 100000){// || an == 1000){
+ toProcess.push_back(std::make_pair(n2t[i1],n2t[i2]));
Msg::Debug("Forbidden quad found in blossom quadrangulation, optimization will be required");
}
else{
@@ -1579,41 +1668,49 @@ static int _recombineIntoQuads(GFace *gf, int recur_level, bool cubicGraph = 1)
return success;
}
-void recombineIntoQuads(GFace *gf)
+void recombineIntoQuads(GFace *gf,
+ bool topologicalOpti,
+ bool nodeRepositioning)
{
double t1 = Cpu();
//gf->model()->writeMSH("before.msh");
- removeFourTrianglesNodes(gf, false);
+ if (topologicalOpti)removeFourTrianglesNodes(gf, false);
int success = _recombineIntoQuads(gf,0);
//gf->model()->writeMSH("raw.msh");
- for(int i = 0; i < CTX::instance()->mesh.nbSmoothing; i++)
- laplaceSmoothing(gf);
+ if (nodeRepositioning)
+ for(int i = 0; i < CTX::instance()->mesh.nbSmoothing; i++)
+ laplaceSmoothing(gf);
- //blossom-quad algo
+ //blossom-quad algo
if(success && CTX::instance()->mesh.algoRecombine == 1){
- //gf->model()->writeMSH("smoothed.msh");
- int COUNT = 0;
- //char NAME[256];
- while(1){
- int x = removeTwoQuadsNodes(gf);
- //if(x){ sprintf(NAME,"iter%dTQ.msh",COUNT++); gf->model()->writeMSH(NAME);}
- int y= removeDiamonds(gf);
- //if(y){ sprintf(NAME,"iter%dD.msh",COUNT++); gf->model()->writeMSH(NAME); }
- int z = _edgeSwapQuadsForBetterQuality ( gf );
+ if (topologicalOpti){
+ //gf->model()->writeMSH("smoothed.msh");
+ int COUNT = 0;
+ //char NAME[256];
+ while(1){
+ int x = removeTwoQuadsNodes(gf);
+ //if(x){ sprintf(NAME,"iter%dTQ.msh",COUNT++); gf->model()->writeMSH(NAME);}
+ int y= removeDiamonds(gf);
+ //if(y){ sprintf(NAME,"iter%dD.msh",COUNT++); gf->model()->writeMSH(NAME); }
+ laplaceSmoothing(gf);
+ int z = 0 ;//edgeSwapQuadsForBetterQuality ( gf );
+ //if(z){ sprintf(NAME,"iter%dS.msh",COUNT++); gf->model()->writeMSH(NAME); }
+ if (!(x+y+z))break;
+ }
+ edgeSwapQuadsForBetterQuality ( gf );
//if(z){ sprintf(NAME,"iter%dS.msh",COUNT++); gf->model()->writeMSH(NAME); }
- if (!(x+y+z))break;
- }
- for(int i = 0; i < CTX::instance()->mesh.nbSmoothing; i++){
- laplaceSmoothing(gf);
+ for(int i = 0; i < CTX::instance()->mesh.nbSmoothing; i++){
+ laplaceSmoothing(gf);
+ }
}
double t2 = Cpu();
Msg::Info("Blossom recombination algorithm completed (%g s)", t2 - t1);
return;
}
- //simple recombination algo
+ //simple recombination algo
_recombineIntoQuads(gf,0);
for(int i = 0; i < CTX::instance()->mesh.nbSmoothing; i++)
laplaceSmoothing(gf);
@@ -1626,3 +1723,37 @@ void recombineIntoQuads(GFace *gf)
Msg::Info("Simple recombination algorithm completed (%g s)", t2 - t1);
}
+
+void quadsToTriangles(GFace *gf){
+ for (int i=0;i<gf->quadrangles.size();i++){
+ MQuadrangle *q = gf->quadrangles[i];
+ MTriangle *t1 = new MTriangle (q->getVertex(0),q->getVertex(1),q->getVertex(2));
+ MTriangle *t2 = new MTriangle (q->getVertex(2),q->getVertex(3),q->getVertex(0));
+ gf->triangles.push_back(t1);
+ gf->triangles.push_back(t2);
+ delete q;
+ }
+ gf->quadrangles.clear();
+}
+
+void recombineIntoQuadsIterative(GFace *gf){
+
+ recombineIntoQuads(gf);
+ // return;
+ int COUNT = 0;
+ while (1){
+ quadsToTriangles(gf);
+ {char NAME[245];sprintf(NAME,"iterT%d.msh",COUNT); gf->model()->writeMSH(NAME);}
+ std::set<MTri3*,compareTri3Ptr> AllTris;
+ std::vector<double> vSizes, vSizesBGM, Us, Vs;
+ std::vector<SMetric3> vMetricsBGM;
+ // buildMeshGenerationDataStructures
+ // (gf, AllTris, vSizes, vSizesBGM, vMetricsBGM, Us, Vs);
+ // int nbSwaps = edgeSwapPass(gf, AllTris, SWCR_DEL, Us, Vs, vSizes, vSizesBGM);
+ // transferDataStructure(gf, AllTris, Us, Vs);
+ {char NAME[245];sprintf(NAME,"iterTD%d.msh",COUNT); gf->model()->writeMSH(NAME);}
+ recombineIntoQuads(gf,false,true);
+ {char NAME[245];sprintf(NAME,"iter%d.msh",COUNT++); gf->model()->writeMSH(NAME);}
+ if (COUNT == 5)break;
+ }
+}
diff --git a/Mesh/meshGFaceOptimize.h b/Mesh/meshGFaceOptimize.h
index bda526530a..327d1d49ec 100644
--- a/Mesh/meshGFaceOptimize.h
+++ b/Mesh/meshGFaceOptimize.h
@@ -30,7 +30,23 @@ 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 buildVertexToElement(std::vector<T*> &eles, v2t_cont &adj){
+ for(unsigned int i = 0; i < eles.size(); i++){
+ T *t = eles[i];
+ for(int j = 0; j < t->getNumVertices(); j++){
+ MVertex *v = t->getVertex(j);
+ v2t_cont :: iterator it = adj.find(v);
+ if(it == adj.end()){
+ std::vector<MElement*> one;
+ one.push_back(t);
+ adj[v] = one;
+ }
+ else{
+ it->second.push_back(t);
+ }
+ }
+ }
+}
template <class T> void buildEdgeToElement(std::vector<T*> &eles, e2t_cont &adj);
void buildVertexToTriangle(std::vector<MTriangle*> &, v2t_cont &adj);
@@ -73,7 +89,10 @@ void buildMeshGenerationDataStructures(GFace *gf,
std::vector<double> &Vs);
void transferDataStructure(GFace *gf, std::set<MTri3*, compareTri3Ptr> &AllTris,
std::vector<double> &Us, std::vector<double> &Vs);
-void recombineIntoQuads(GFace *gf);
+void recombineIntoQuads(GFace *gf,
+ bool topologicalOpti = true,
+ bool nodeRepositioning = true);
+void recombineIntoQuadsIterative(GFace *gf);
struct swapquad{
int v[4];
diff --git a/Mesh/meshGFaceQuadrilateralize.cpp b/Mesh/meshGFaceQuadrilateralize.cpp
index 27865202e7..ce2bd31e7f 100644
--- a/Mesh/meshGFaceQuadrilateralize.cpp
+++ b/Mesh/meshGFaceQuadrilateralize.cpp
@@ -292,7 +292,8 @@ bool edgeFront::formQuad(BDS_Edge *e, BDS_Edge *left, BDS_Edge *right)
// printf("recover\n");
// top = m->recover_edge_fast(pleft,pright);
// if(!top)
- top = m->recover_edge(pleft->iD, pright->iD);
+ bool _fatallyFailed;
+ top = m->recover_edge(pleft->iD, pright->iD,_fatallyFailed);
// printf("recover done %p\n",top);
if(!top) return false;
}
diff --git a/Mesh/meshRefine.cpp b/Mesh/meshRefine.cpp
index 496834387a..4cefa59b8e 100644
--- a/Mesh/meshRefine.cpp
+++ b/Mesh/meshRefine.cpp
@@ -377,8 +377,10 @@ 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)
- Subdivide(*it, splitIntoQuads, splitIntoHexas, faceVertices);
+ for(GModel::fiter it = m->firstFace(); it != m->lastFace(); ++it){
+ bool splitIntoQuadsForThisFace = true;
+ Subdivide(*it, splitIntoQuadsForThisFace, splitIntoHexas, faceVertices);
+ }
for(GModel::riter it = m->firstRegion(); it != m->lastRegion(); ++it)
Subdivide(*it, splitIntoHexas, faceVertices);
diff --git a/benchmarks/boolean/sphere.lua b/benchmarks/boolean/sphere.lua
index 78ace91695..6477c57e48 100644
--- a/benchmarks/boolean/sphere.lua
+++ b/benchmarks/boolean/sphere.lua
@@ -1,6 +1,6 @@
options = gmshOptions()
options:initOptions()
-options:numberSet('Mesh', 0, 'CharacteristicLengthFactor', 0.9grep -IIr)
+options:numberSet('Mesh', 0, 'CharacteristicLengthFactor', 0.9);
R = 1.0;
myModel = GModel();
@@ -8,4 +8,4 @@ myModel:addSphere(0,0,0,R);
myModel:setAsCurrent();
-myModel:mesh(2);
+--myModel:mesh(2);
--
GitLab