From 2a34f6f5f251fb86e07858c40b0f6ed02397c402 Mon Sep 17 00:00:00 2001
From: Jean-Francois Remacle <jean-francois.remacle@uclouvain.be>
Date: Tue, 21 Sep 2010 14:25:51 +0000
Subject: [PATCH] meshes of spheres should work well !!!
---
Geo/OCCFace.cpp | 9 ++++++
Mesh/BDS.cpp | 22 ++++++++++++--
Mesh/meshGFaceBDS.cpp | 56 ++++++++++++++++++++++++-----------
Mesh/meshGFaceOptimize.cpp | 21 ++++++++++---
benchmarks/boolean/sphere.lua | 13 +++-----
benchmarks/step/tank.geo | 1 +
6 files changed, 90 insertions(+), 32 deletions(-)
diff --git a/Geo/OCCFace.cpp b/Geo/OCCFace.cpp
index 55dbff4305..f9619ef445 100644
--- a/Geo/OCCFace.cpp
+++ b/Geo/OCCFace.cpp
@@ -96,6 +96,15 @@ void OCCFace::setup()
occface = BRep_Tool::Surface(s);
// specific parametrization for a sphere.
_isSphere = isSphere(_radius,_center);
+ if (_isSphere){
+ for (std::list<GEdge*>::iterator it = l_edges.begin();
+ it != l_edges.end(); ++it){
+ GEdge *ge = *it;
+ if (ge->isSeam(this)){
+
+ }
+ }
+ }
}
Range<double> OCCFace::parBounds(int i) const
diff --git a/Mesh/BDS.cpp b/Mesh/BDS.cpp
index ec602f3363..1c98232a95 100644
--- a/Mesh/BDS.cpp
+++ b/Mesh/BDS.cpp
@@ -1156,6 +1156,7 @@ bool BDS_Mesh::collapse_edge_parametric(BDS_Edge *e, BDS_Point *p)
bool BDS_Mesh::smooth_point_centroid(BDS_Point *p, GFace *gf, bool test_quality)
{
+
if(!p->config_modified) return false;
if(p->g && p->g->classif_degree <= 1) return false;
if(p->g && p->g->classif_tag < 0) {
@@ -1168,6 +1169,10 @@ bool BDS_Mesh::smooth_point_centroid(BDS_Point *p, GFace *gf, bool test_quality)
eit++;
}
+ /* TEST */
+ double R; SPoint3 c; bool isSphere = gf->isSphere(R,c);
+ double XX=0,YY=0,ZZ=0;
+
double U = 0;
double V = 0;
double LC = 0;
@@ -1187,14 +1192,27 @@ bool BDS_Mesh::smooth_point_centroid(BDS_Point *p, GFace *gf, bool test_quality)
sTot += fact;
U += n->u * fact;
V += n->v * fact;
+ XX += n->X;
+ YY += n->Y;
+ ZZ += n->Z;
LC += n->lc() * fact;
++ited;
}
U /= (sTot);
V /= (sTot);
LC /= (sTot);
-
- GPoint gp = gf->point(U * scalingU, V * scalingV);
+ XX/= (sTot);
+ YY/= (sTot);
+ ZZ/= (sTot);
+
+ GPoint gp;double uv[2];
+ if (isSphere){
+ gp = gf->closestPoint(SPoint3(XX, YY, ZZ), uv);
+ U = gp.u();
+ V = gp.v();
+ }
+ else
+ gp = gf->point(U * scalingU, V * scalingV);
if (!gp.succeeded()){
return false;
diff --git a/Mesh/meshGFaceBDS.cpp b/Mesh/meshGFaceBDS.cpp
index b4e0ed943e..5a63eb5125 100644
--- a/Mesh/meshGFaceBDS.cpp
+++ b/Mesh/meshGFaceBDS.cpp
@@ -451,35 +451,52 @@ void delaunayizeBDS(GFace *gf, BDS_Mesh &m, int &nb_swap)
}
// A test for spheres
-static void midpointsphere(GFace *gf, double u1, double v1, double u2,
- double v2, double &u12, double &v12, double r)
+static void midpointsphere(GFace *gf,
+ double u1,
+ double v1,
+ double u2,
+ double v2,
+ double &u12,
+ double &v12,
+ SPoint3 ¢er,
+ double r)
{
GPoint p1 = gf->point(u1, v1);
GPoint p2 = gf->point(u2, v2);
+
+ SVector3 DIR ((p1.x()+p2.x())/2.0 - center.x(),
+ (p1.y()+p2.y())/2.0 - center.y(),
+ (p1.z()+p2.z())/2.0 - center.z());
+ DIR.normalize();
+
+ double X,Y,Z;
+
double guess [2] = {0.5 * (u1 + u2), 0.5 * (v1 + v2)};
u12 = guess[0];
v12 = guess[1];
-
- double d = sqrt((p1.x() - p2.x()) * (p1.x() - p2.x()) +
- (p1.y() - p2.y()) * (p1.y() - p2.y()) +
- (p1.z() - p2.z()) * (p1.z() - p2.z()));
-
- if (d > r/3){
+ if ( (v1 > 0.7*M_PI/2 && v2 > 0.7*M_PI/2) ||
+ (v1 < -0.7*M_PI/2 && v2 < -0.7*M_PI/2) ){
+ // printf("coucou\n");
+ X = center.x() + DIR.x() * r;
+ Y = center.y() + DIR.y() * r;
+ Z = center.z() + DIR.z() * r;
+ }
+ else{
return;
+ X = center.x() - DIR.x() * r;
+ Y = center.y() - DIR.y() * r;
+ Z = center.z() - DIR.z() * r;
}
- const double X = 0.5 * (p1.x() + p2.x());
- const double Y = 0.5 * (p1.y() + p2.y());
- const double Z = 0.5 * (p1.z() + p2.z());
-
+
GPoint proj = gf->closestPoint(SPoint3(X, Y, Z), guess);
if (proj.succeeded()){
u12 = proj.u();
v12 = proj.v();
}
-
+ // printf("%g %g %g -- %g\n",center.x(),center.y(),center.z(),r);
+ // printf("%g %g -- %g %g -- %g %g -- %g %g\n",
+ // u1, v1, u2, v2, u12, v12, 0.5 * (u1 + u2), 0.5 * (v1 + v2));
return;
- printf("%g %g -- %g %g -- %g %g -- %g %g\n",
- u1, v1, u2, v2, u12, v12, 0.5 * (u1 + u2), 0.5 * (v1 + v2));
}
void splitEdgePass(GFace *gf, BDS_Mesh &m, double MAXE_, int &nb_split)
@@ -499,6 +516,10 @@ void splitEdgePass(GFace *gf, BDS_Mesh &m, double MAXE_, int &nb_split)
std::sort(edges.begin(), edges.end());
+ double RADIUS;
+ SPoint3 CENTER;
+ bool isSphere = gf->isSphere(RADIUS,CENTER);
+
for (unsigned int i = 0; i < edges.size(); ++i){
BDS_Edge *e = edges[i].second;
if (!e->deleted){
@@ -508,9 +529,9 @@ void splitEdgePass(GFace *gf, BDS_Mesh &m, double MAXE_, int &nb_split)
BDS_Point *mid ;
double U, V;
- if (0 && gf->geomType() == GEntity::Sphere){
+ if (0 && isSphere){
midpointsphere(gf,e->p1->u,e->p1->v,e->p2->u,e->p2->v,U,V,
- gf-> getSurfaceParams().radius);
+ CENTER,RADIUS);
}
else{
U = coord * e->p1->u + (1 - coord) * e->p2->u;
@@ -603,6 +624,7 @@ void collapseEdgePassUnSorted(GFace *gf, BDS_Mesh &m, double MINE_, int MAXNP,
void smoothVertexPass(GFace *gf, BDS_Mesh &m, int &nb_smooth, bool q)
{
+ // return;
std::set<BDS_Point*,PointLessThan>::iterator itp = m.points.begin();
while(itp != m.points.end()){
if(m.smooth_point_centroid(*itp, gf,q))
diff --git a/Mesh/meshGFaceOptimize.cpp b/Mesh/meshGFaceOptimize.cpp
index e2a7d50156..c93707ec0c 100644
--- a/Mesh/meshGFaceOptimize.cpp
+++ b/Mesh/meshGFaceOptimize.cpp
@@ -787,6 +787,9 @@ void laplaceSmoothing(GFace *gf)
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()){
@@ -799,27 +802,37 @@ void laplaceSmoothing(GFace *gf)
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, cu / fact);
- ver->setParameter(1, cv / fact);
- GPoint pt = gf->point(SPoint2(cu / fact, cv / fact));
+ 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();
diff --git a/benchmarks/boolean/sphere.lua b/benchmarks/boolean/sphere.lua
index f179218ca6..a826c4c751 100644
--- a/benchmarks/boolean/sphere.lua
+++ b/benchmarks/boolean/sphere.lua
@@ -1,10 +1,5 @@
+R = 1.0;
+myModel = GModel();
+myModel:addSphere(0,0,0,R);
-options = gmshOptions()
-options:numberSet('Mesh', 0, 'CharacteristicLengthFactor', 1.5)
-
-myTool = GModel();
-myTool:addSphere(0.0,0.0,0.0,1);
-
-myTool:mesh(1)
-myTool:mesh(2)
-myTool:save('sphere.msh')
+myModel:setAsCurrent();
diff --git a/benchmarks/step/tank.geo b/benchmarks/step/tank.geo
index 123aefd78d..513a3baaac 100644
--- a/benchmarks/step/tank.geo
+++ b/benchmarks/step/tank.geo
@@ -40,3 +40,4 @@ Transfinite Surface {14} = {21,25,26,18};
Mesh.Smoothing = 10;
Recombine Surface{1:52};
+Mesh.SubdivisionAlgorithm=3;
\ No newline at end of file
--
GitLab