From 656278c91c71e240d7ffd4fdc9f4a1f60703d8e6 Mon Sep 17 00:00:00 2001
From: Jean-Francois Remacle <jean-francois.remacle@uclouvain.be>
Date: Tue, 26 Feb 2013 13:30:50 +0000
Subject: [PATCH] spallart almaras model gets better, bug fixed in delaunay
surface meshing with embedded edges bug fixed in cross filed computation for
non conformal mappings
---
CMakeLists.txt | 6 +-
Geo/SVector3.h | 8 +++
Geo/intersectCurveSurface.cpp | 2 +
Mesh/BackgroundMesh.cpp | 103 ++++++++++++----------------
Mesh/meshGFace.cpp | 27 --------
Mesh/meshGFaceDelaunayInsertion.cpp | 39 +++++------
Mesh/meshGFaceDelaunayInsertion.h | 1 +
Mesh/meshGFaceOptimize.cpp | 24 ++++---
Mesh/surfaceFiller.cpp | 56 +++++++++++----
benchmarks/hex/twoHoles.geo | 2 +-
10 files changed, 136 insertions(+), 132 deletions(-)
diff --git a/CMakeLists.txt b/CMakeLists.txt
index 99b4ee3b68..cbd192b676 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -563,11 +563,13 @@ endif(ENABLE_MPI)
if(ENABLE_POPPLER)
find_library(POPPLER_LIB poppler)
find_library(POPPLER_CPP_LIB poppler-cpp)
- if(POPPLER_LIB)
+ find_path(POPPLER_INC "poppler/cpp/poppler-document.h" PATH_SUFFIXES src include)
+ if(POPPLER_LIB AND POPPLER_INC)
set_config_option(HAVE_POPPLER "Poppler")
list(APPEND EXTERNAL_LIBRARIES ${POPPLER_LIB})
list(APPEND EXTERNAL_LIBRARIES ${POPPLER_CPP_LIB})
- endif(POPPLER_LIB)
+ list(APPEND EXTERNAL_INCLUDES ${POPPLER_INC})
+ endif(POPPLER_LIB AND POPPLER_INC)
endif(ENABLE_POPPLER)
if(HAVE_MESH OR HAVE_SOLVER)
diff --git a/Geo/SVector3.h b/Geo/SVector3.h
index e42e2f9cc7..f58a069cf0 100644
--- a/Geo/SVector3.h
+++ b/Geo/SVector3.h
@@ -89,6 +89,13 @@ inline SVector3 crossprod(const SVector3 &a, const SVector3 &b)
-(a.x() * b.z() - b.x() * a.z()),
a.x() * b.y() - b.x() * a.y()); }
+inline double angle (const SVector3 &a, const SVector3 &b){
+ double cosTheta = dot(a,b);
+ double sinTheta = norm(crossprod(a,b));
+ return atan2 (sinTheta,cosTheta);
+}
+
+
inline SVector3 operator*(double m,const SVector3 &v)
{ return SVector3(v[0] * m, v[1] * m, v[2] * m); }
@@ -104,6 +111,7 @@ inline SVector3 operator+(const SVector3 &a,const SVector3 &b)
inline SVector3 operator-(const SVector3 &a,const SVector3 &b)
{ return SVector3(a[0] - b[0], a[1] - b[1], a[2] - b[2]); }
+
inline void buildOrthoBasis_naive(SVector3 &dir, SVector3 &dir1, SVector3 &dir2)
{
dir.normalize();
diff --git a/Geo/intersectCurveSurface.cpp b/Geo/intersectCurveSurface.cpp
index 9e51a3790d..f94b276f44 100644
--- a/Geo/intersectCurveSurface.cpp
+++ b/Geo/intersectCurveSurface.cpp
@@ -21,6 +21,7 @@ struct intersectCurveSurfaceData
uvt(2) = newPoint[2];
fullVector<double> res(3);
_kaboom(uvt,res,this);
+ // printf("start with %12.5E\n",res.norm());
if (res.norm() < epsilon)return true;
if(newton_fd(_kaboom, uvt, this)){
@@ -41,6 +42,7 @@ void _kaboom(fullVector<double> &uvt,
res(0) = s.x() - c.x();
res(1) = s.y() - c.y();
res(2) = s.z() - c.z();
+ // printf("%g %g %g vs %g %g %g \n",s.x(),s.y(),s.z(),c.x(),c.y(),c.z());
}
int intersectCurveSurface (curveFunctor &c, surfaceFunctor & s, double uvt[3], double epsilon){
diff --git a/Mesh/BackgroundMesh.cpp b/Mesh/BackgroundMesh.cpp
index 013a45b541..3dd3273077 100644
--- a/Mesh/BackgroundMesh.cpp
+++ b/Mesh/BackgroundMesh.cpp
@@ -623,8 +623,9 @@ static void propagateValuesOnFace(GFace *_gf,
}
// Solve
- if (myAssembler.sizeOfR())
+ if (myAssembler.sizeOfR()){
_lsys->systemSolve();
+ }
// save solution
for (std::set<MVertex*>::iterator it = vs.begin(); it != vs.end(); ++it){
@@ -653,18 +654,19 @@ void backgroundMesh::propagate1dMesh(GFace *_gf)
for(unsigned int i = 0; i < (*it)->lines.size(); i++ ){
MVertex *v1 = (*it)->lines[i]->getVertex(0);
MVertex *v2 = (*it)->lines[i]->getVertex(1);
- // printf("%g %g %g\n",v1->x(),v1->y(),v1->z());
- double d = sqrt((v1->x() - v2->x()) * (v1->x() - v2->x()) +
- (v1->y() - v2->y()) * (v1->y() - v2->y()) +
- (v1->z() - v2->z()) * (v1->z() -v2->z()));
- for (int k=0;k<2;k++){
- MVertex *v = (*it)->lines[i]->getVertex(k);
- std::map<MVertex*, double>::iterator itv = sizes.find(v);
- if (itv == sizes.end())
- sizes[v] = log(d);
- else
- itv->second = 0.5 * (itv->second + log(d));
- }
+ if (v1 != v2){
+ double d = sqrt((v1->x() - v2->x()) * (v1->x() - v2->x()) +
+ (v1->y() - v2->y()) * (v1->y() - v2->y()) +
+ (v1->z() - v2->z()) * (v1->z() -v2->z()));
+ for (int k=0;k<2;k++){
+ MVertex *v = (*it)->lines[i]->getVertex(k);
+ std::map<MVertex*, double>::iterator itv = sizes.find(v);
+ if (itv == sizes.end())
+ sizes[v] = log(d);
+ else
+ itv->second = 0.5 * (itv->second + log(d));
+ }
+ }
}
}
}
@@ -711,19 +713,26 @@ void backgroundMesh::propagateCrossFieldByDistance(GFace *_gf)
v[1] = (*it)->lines[i]->getVertex(1);
SPoint2 p1,p2;
reparamMeshEdgeOnFace(v[0],v[1],_gf,p1,p2);
- double angle = atan2 ( p1.y()-p2.y() , p1.x()-p2.x() );
- crossField2d::normalizeAngle (angle);
+ /* a correct way of computing angles */
+ Pair<SVector3, SVector3> der = _gf->firstDer((p1+p2)*.5);
+ SVector3 t1 = der.first();
+ SVector3 t2 (v[1]->x()-v[0]->x(),v[1]->y()-v[0]->y(),v[1]->z()-v[0]->z());
+ t1.normalize();
+ t2.normalize();
+ double _angle = angle (t1,t2);
+ // double angle = atan2 ( p1.y()-p2.y() , p1.x()-p2.x() );
+ crossField2d::normalizeAngle (_angle);
for (int i=0;i<2;i++){
std::map<MVertex*,double>::iterator itc = _cosines4.find(v[i]);
std::map<MVertex*,double>::iterator its = _sines4.find(v[i]);
if (itc != _cosines4.end()){
- itc->second = 0.5*(itc->second + cos(4*angle));
- its->second = 0.5*(its->second + sin(4*angle));
+ itc->second = 0.5*(itc->second + cos(4*_angle));
+ its->second = 0.5*(its->second + sin(4*_angle));
}
else {
_param[v[i]] = (i==0) ? p1 : p2;
- _cosines4[v[i]] = cos(4*angle);
- _sines4[v[i]] = sin(4*angle);
+ _cosines4[v[i]] = cos(4*_angle);
+ _sines4[v[i]] = sin(4*_angle);
}
}
}
@@ -756,6 +765,7 @@ void backgroundMesh::propagateCrossFieldByDistance(GFace *_gf)
void backgroundMesh::propagatecrossField(GFace *_gf)
{
+
std::map<MVertex*,double> _cosines4,_sines4;
std::list<GEdge*> e;
@@ -771,51 +781,27 @@ void backgroundMesh::propagatecrossField(GFace *_gf)
v[1] = (*it)->lines[i]->getVertex(1);
SPoint2 p1,p2;
reparamMeshEdgeOnFace(v[0],v[1],_gf,p1,p2);
- double angle = atan2 ( p1.y()-p2.y() , p1.x()-p2.x() );
+ Pair<SVector3, SVector3> der = _gf->firstDer((p1+p2)*.5);
+ SVector3 t1 = der.first();
+ SVector3 t2 (v[1]->x()-v[0]->x(),v[1]->y()-v[0]->y(),v[1]->z()-v[0]->z());
+ t1.normalize();
+ t2.normalize();
+ double _angle = angle (t1,t2);
+ // printf("angle = %12.5E\n",_angle);
+ // angle_ = atan2 ( p1.y()-p2.y() , p1.x()-p2.x() );
//double angle = atan2 ( v[0]->y()-v[1]->y() , v[0]->x()- v[1]->x() );
//double angle = atan2 ( v0->y()-v1->y() , v0->x()- v1->x() );
- crossField2d::normalizeAngle (angle);
+ crossField2d::normalizeAngle (_angle);
for (int i=0;i<2;i++){
std::map<MVertex*,double>::iterator itc = _cosines4.find(v[i]);
std::map<MVertex*,double>::iterator its = _sines4.find(v[i]);
if (itc != _cosines4.end()){
- itc->second = 0.5*(itc->second + cos(4*angle));
- its->second = 0.5*(its->second + sin(4*angle));
+ itc->second = 0.5*(itc->second + cos(4*_angle));
+ its->second = 0.5*(its->second + sin(4*_angle));
}
else {
- _cosines4[v[i]] = cos(4*angle);
- _sines4[v[i]] = sin(4*angle);
- }
- }
- }
- }
- }
-
- // force smooth transition
- const int nbSmooth = 0;
- const double threshold_angle = 2. * M_PI/180.;
- for (int SMOOTH_ITER = 0 ; SMOOTH_ITER < nbSmooth ; SMOOTH_ITER++){
- it = e.begin();
- for( ; it != e.end(); ++it ){
- if (!(*it)->isSeam(_gf)){
- for(unsigned int i = 0; i < (*it)->lines.size(); i++ ){
- MVertex *v[2];
- v[0] = (*it)->lines[i]->getVertex(0);
- v[1] = (*it)->lines[i]->getVertex(1);
- double cos40 = _cosines4[v[0]];
- double cos41 = _cosines4[v[1]];
- double sin40 = _sines4[v[0]];
- double sin41 = _sines4[v[1]];
- double angle0 = atan2 (sin40,cos40)/4.;
- double angle1 = atan2 (sin41,cos41)/4.;
- if (fabs(angle0 - angle1) > threshold_angle ){
- double angle0_new = angle0 - (angle0-angle1) * 0.1;
- double angle1_new = angle1 + (angle0-angle1) * 0.1;
- // printf("%g %g -- %g %g\n",angle0,angle1,angle0_new,angle1_new);
- _cosines4[v[0]] = cos(4*angle0_new);
- _sines4[v[0]] = sin(4*angle0_new);
- _cosines4[v[1]] = cos(4*angle1_new);
- _sines4[v[1]] = sin(4*angle1_new);
+ _cosines4[v[i]] = cos(4*_angle);
+ _sines4[v[i]] = sin(4*_angle);
}
}
}
@@ -1013,7 +999,7 @@ void backgroundMesh::print(const std::string &filename, GFace *gf,
*/
}
else {
- /*
+
GPoint p1 = gf->point(SPoint2(v1->x(),v1->y()));
GPoint p2 = gf->point(SPoint2(v2->x(),v2->y()));
GPoint p3 = gf->point(SPoint2(v3->x(),v3->y()));
@@ -1021,12 +1007,13 @@ void backgroundMesh::print(const std::string &filename, GFace *gf,
p1.x(),p1.y(),p1.z(),
p2.x(),p2.y(),p2.z(),
p3.x(),p3.y(),p3.z(),itv1->second,itv2->second,itv3->second);
- */
+ /*
fprintf(f,"ST(%g,%g,%g,%g,%g,%g,%g,%g,%g) {%g,%g,%g};\n",
v1->x(),v1->y(),v1->z(),
v2->x(),v2->y(),v2->z(),
v3->x(),v3->y(),v3->z(),
itv1->second,itv2->second,itv3->second);
+ */
}
}
fprintf(f,"};\n");
diff --git a/Mesh/meshGFace.cpp b/Mesh/meshGFace.cpp
index fd35db2185..ce37386e01 100644
--- a/Mesh/meshGFace.cpp
+++ b/Mesh/meshGFace.cpp
@@ -101,45 +101,18 @@ static void copyMesh(GFace *source, GFace *target)
}
}
- SPoint2 param_source[2], param_target[2];
SVector3 DX;
int count = 0;
for (std::map<MVertex*, MVertex*>::iterator it = vs2vt.begin(); it != vs2vt.end() ; ++it){
MVertex *vs = it->first;
MVertex *vt = it->second;
- reparamMeshVertexOnFace(vs, source, param_source[count]);
- reparamMeshVertexOnFace(vt, target, param_target[count++]);
DX = SVector3(vt->x() - vs->x(), vt->y() - vs->y(), vt->z() - vs->z());
if (count == 2) break;
}
- // double t1u = param_target[0].x(), t1v = param_target[0].y();
- // double t2u = param_target[1].x(), t2v = param_target[1].y();
- // double s1u = param_source[0].x(), s1v = param_source[0].y();
- // double s2u = param_source[1].x(), s2v = param_source[1].y();
-
- // SVector3 _a(s2u - s1u, s2v - s1v, 0);
- // SVector3 _b(t2u - t1u, t2v - t1v, 0);
-
- // SVector3 _c = crossprod(_a, _b);
- // double sinA = _c.z();
- // double cosA = dot(_a, _b);
- // const double theta = atan2(sinA, cosA);
- // const double c = cos(theta);
- // const double s = sin(theta);
-
for(unsigned int i = 0; i < source->mesh_vertices.size(); i++){
MVertex *vs = source->mesh_vertices[i];
- double u, v;
- vs->getParameter(0, u);
- vs->getParameter(1, v);
- // apply transformation
- // const double U = c * (u - s1u) + s * (v - s1v) + t1u;
- // const double V = -s * (u - s1u) + c * (v - s1v) + t1v;
SPoint3 tp (vs->x() + DX.x(),vs->y() + DX.y(),vs->z() + DX.z());
- // const double initialGuess[2] = {U,V};
- // FIXME !!!
- // assume a translation for now !!!
SPoint2 XXX = target->parFromPoint(tp);
GPoint gp = target->point(XXX);
diff --git a/Mesh/meshGFaceDelaunayInsertion.cpp b/Mesh/meshGFaceDelaunayInsertion.cpp
index 6adf207b67..b1d210bf92 100644
--- a/Mesh/meshGFaceDelaunayInsertion.cpp
+++ b/Mesh/meshGFaceDelaunayInsertion.cpp
@@ -490,15 +490,18 @@ void recurFindCavityAniso(GFace *gf,
for (int i = 0; i < 3; i++){
MTri3 *neigh = t->getNeigh(i) ;
- if (!neigh)
- shell.push_back(edgeXface(t, i));
+ edgeXface exf (t, i);
+ // take care of untouchable internal edges
+ std::set<MEdge,Less_Edge>::iterator it = data.internalEdges.find(MEdge(exf.v[0],exf.v[1]));
+ if (!neigh || it != data.internalEdges.end())
+ shell.push_back(exf);
else if (!neigh->isDeleted()){
// int circ = inCircumCircleAniso(gf, neigh->tri(), param, metric, data);
int circ = inCircumCircleAniso(gf, neigh->tri(), param, metric, data);
if (circ)
recurFindCavityAniso(gf, shell, cavity,metric, param, neigh, data);
else
- shell.push_back(edgeXface(t, i));
+ shell.push_back(exf);
}
}
}
@@ -791,30 +794,26 @@ static bool insertAPoint(GFace *gf, std::set<MTri3*,compareTri3Ptr>::iterator it
double uv[2];
// if the point is able to break the bad triangle "worst"
- if (1){
- if (inCircumCircleAniso(gf, worst->tri(), center, metric, data)){
- // double t1 = Cpu();
- recurFindCavityAniso(gf, shell, cavity, metric, center, worst, data);
- // __DT1 += Cpu() - t1 ;
- for (std::list<MTri3*>::iterator itc = cavity.begin(); itc != cavity.end(); ++itc){
- if (invMapUV((*itc)->tri(), center, data, uv, 1.e-8)) {
- ptin = *itc;
- break;
- }
+ if (inCircumCircleAniso(gf, worst->tri(), center, metric, data)){
+ // double t1 = Cpu();
+ recurFindCavityAniso(gf, shell, cavity, metric, center, worst, data);
+ // __DT1 += Cpu() - t1 ;
+ for (std::list<MTri3*>::iterator itc = cavity.begin(); itc != cavity.end(); ++itc){
+ if (invMapUV((*itc)->tri(), center, data, uv, 1.e-8)) {
+ ptin = *itc;
+ break;
}
}
- // else look for it
- else {
- // printf("cocuou\n");
- ptin = search4Triangle (worst, center, data, AllTris,uv, oneNewTriangle ? true : false);
+ }
+ // else look for it
+ else {
+ // printf("cocuou\n");
+ ptin = search4Triangle (worst, center, data, AllTris,uv, oneNewTriangle ? true : false);
if (ptin) {
recurFindCavityAniso(gf, shell, cavity, metric, center, ptin, data);
}
- }
}
- // ptin = search4Triangle (worst, center, Us, Vs, AllTris,uv, oneNewTriangle ? true : false);
-
if (ptin) {
// we use here local coordinates as real coordinates
// x,y and z will be computed hereafter
diff --git a/Mesh/meshGFaceDelaunayInsertion.h b/Mesh/meshGFaceDelaunayInsertion.h
index 00ab44a4c4..96590d8c08 100644
--- a/Mesh/meshGFaceDelaunayInsertion.h
+++ b/Mesh/meshGFaceDelaunayInsertion.h
@@ -26,6 +26,7 @@ struct bidimMeshData
std::vector<SMetric3> vMetricsBGM;
std::map<MVertex* , MVertex*>* equivalence;
std::map<MVertex*, SPoint2> * parametricCoordinates;
+ std::set<MEdge,Less_Edge> internalEdges; // embedded edges
inline void addVertex (MVertex* mv, double u, double v, double size, double sizeBGM){
int index = Us.size();
if (mv->onWhat()->dim() == 2)mv->setIndex(index);
diff --git a/Mesh/meshGFaceOptimize.cpp b/Mesh/meshGFaceOptimize.cpp
index e2653544fa..c27522a364 100644
--- a/Mesh/meshGFaceOptimize.cpp
+++ b/Mesh/meshGFaceOptimize.cpp
@@ -163,21 +163,25 @@ void buildMeshGenerationDataStructures(GFace *gf,
}
}
+ // take good care of embedded edges
+ {
+ std::list<GEdge*> embedded_edges = gf->embeddedEdges();
+ std::list<GEdge*>::iterator ite = embedded_edges.begin();
+ while(ite != embedded_edges.end()){
+ if(!(*ite)->isMeshDegenerated()){
+ for (int i=0;i<(*ite)->lines.size();i++)
+ data.internalEdges.insert (MEdge((*ite)->lines[i]->getVertex(0),
+ (*ite)->lines[i]->getVertex(1)));
+ }
+ ++ite;
+ }
+ }
+
// int NUM = 0;
for(std::map<MVertex*, double>::iterator it = vSizesMap.begin();
it != vSizesMap.end(); ++it){
- // FIXME: this vertex-stored indexing makes the optimization
- // routines not thread-safe (we cannot concurrently optimize two
- // surfaces that share an edge)
- // it->first->setIndex(NUM++);
- // OK, I can fix that
- // vSizes.push_back(it->second);
- // vSizesBGM.push_back(it->second);
- // vMetricsBGM.push_back(SMetric3(it->second));
SPoint2 param;
reparamMeshVertexOnFace(it->first, gf, param);
- // Us.push_back(param[0]);
- // Vs.push_back(param[1]);
data.addVertex (it->first, param[0], param[1], it->second, it->second);
}
for(unsigned int i = 0; i < gf->triangles.size(); i++){
diff --git a/Mesh/surfaceFiller.cpp b/Mesh/surfaceFiller.cpp
index 58690ca631..490bdccdbb 100644
--- a/Mesh/surfaceFiller.cpp
+++ b/Mesh/surfaceFiller.cpp
@@ -142,7 +142,7 @@ bool inExclusionZone (SPoint2 &p,
if (!backgroundMesh::current()->inDomain(p.x(),p.y(),0)) return true;
my_wrapper w (p);
- double _min[2] = {p.x()-1.e-3, p.y()-1.e-3},_max[2] = {p.x()+1.e-3,p.y()+1.e-3};
+ double _min[2] = {p.x()-1.e-1, p.y()-1.e-1},_max[2] = {p.x()+1.e-1,p.y()+1.e-1};
rtree.Search(_min,_max,rtree_callback,&w);
return w._tooclose;
@@ -239,11 +239,13 @@ bool compute4neighbors (GFace *gf, // the surface
// and get covariant coordinates a and b
double rhs1[2] = {dot(t1,s1),dot(t1,s2)}, covar1[2];
if (!sys2x2(metric,rhs1,covar1)){
+ Msg::Info("Argh");
covar1[1] = 1.0; covar1[0] = 0.0;
}
double rhs2[2] = {dot(t2,s1),dot(t2,s2)}, covar2[2];
if (!sys2x2(metric,rhs2,covar2)){
- covar2[0] = 1.0; covar2[1] = 0.0;
+ Msg::Info("Argh");
+ covar2[0] = 1.0; covar2[1] = 0.0;
}
// transform the sizes with respect to the metric
@@ -264,7 +266,7 @@ bool compute4neighbors (GFace *gf, // the surface
2*E*covar2[1]*covar2[0]+
N*covar2[1]*covar2[1]);
- // if (l1 == 0.0 || l2 == 0.0) printf("bouuuuuuuuuuuuh %g %g %g %g --- %g %g %g %g %g %g\n",l1,l2,t1.norm(),t2.norm(),s1.x(),s1.y(),s1.z(),s2.x(),s2.y(),s2.z());
+ //if (l1 == 0.0 || l2 == 0.0) printf("bouuuuuuuuuuuuh %g %g %g %g --- %g %g %g %g %g %g\n",l1,l2,t1.norm(),t2.norm(),s1.x(),s1.y(),s1.z(),s2.x(),s2.y(),s2.z());
/* printf("%12.5E %12.5E %12.5E %12.5E %12.5E %12.5E %12.5E %12.5E %g %g %g %g %g %g %g %g %g %g %g\n",
M*covar1[0]*covar1[0]+
@@ -296,18 +298,42 @@ bool compute4neighbors (GFace *gf, // the surface
// We could stop here. Yet, if the metric varies a lot, we can solve
// a nonlinear problem in order to find a better approximation in the real
// surface
+ double ERR[4];
+ for (int i=0;i<4;i++){ //
+ GPoint pp = gf->point(SPoint2(newPoint[i][0], newPoint[i][1]));
+ double D = sqrt ((pp.x() - v_center->x())*(pp.x() - v_center->x()) +
+ (pp.y() - v_center->y())*(pp.y() - v_center->y()) +
+ (pp.z() - v_center->z())*(pp.z() - v_center->z()) );
+ ERR[i] = 100*fabs(D-L)/(D+L);
+ // printf("L = %12.5E D = %12.5E ERR = %12.5E\n",L,D,100*fabs(D-L)/(D+L));
+ }
+
if (goNonLinear){//---------------------------------------------------//
surfaceFunctorGFace ss (gf); //
- SVector3 dirs[4] = {t1*(-1.0),t2*(-1.0),t1,t2}; //
+ SVector3 dirs[4] = {t1*(-1.0),t2*(-1.0),t1*(1.0),t2*(1.0)}; //
for (int i=0;i<4;i++){ //
- double uvt[3] = {newPoint[i][0],newPoint[i][1],0.0}; //
- curveFunctorCircle cf (n,dirs[i],
- SVector3(v_center->x(),v_center->y(),v_center->z()),
- (i%2==1 ?size_param_1:size_param_2)*.5); //
- if (intersectCurveSurface (cf,ss,uvt,size_param_1*1.e-8)){ //
- newPoint[i][0] = uvt[0]; //
- newPoint[i][1] = uvt[1]; //
- } //
+ if (ERR[i] > 12){
+ double uvt[3] = {newPoint[i][0],newPoint[i][1],0.0}; //
+ // printf("Intersecting with circle N = %g %g %g dir = %g %g %g R = %g p = %g %g %g\n",n.x(),n.y(),n.z(),dirs[i].x(),dirs[i].y(),dirs[i].z(),L,v_center->x(),v_center->y(),v_center->z());
+ curveFunctorCircle cf (dirs[i],n,
+ SVector3(v_center->x(),v_center->y(),v_center->z()),
+ L);
+ if (intersectCurveSurface (cf,ss,uvt,size_param_1*1.e-5)){ //
+ GPoint pp = gf->point(SPoint2(uvt[0],uvt[1]));
+ double D = sqrt ((pp.x() - v_center->x())*(pp.x() - v_center->x()) +
+ (pp.y() - v_center->y())*(pp.y() - v_center->y()) +
+ (pp.z() - v_center->z())*(pp.z() - v_center->z()) );
+ double newErr = 100*fabs(D-L)/(D+L);
+ if (newErr < 1){
+ // printf("%12.5E vs %12.5E : %12.5E %12.5E vs %12.5E %12.5E \n",ERR[i],newErr,newPoint[i][0],newPoint[i][1],uvt[0],uvt[1]);
+ newPoint[i][0] = uvt[0]; //
+ newPoint[i][1] = uvt[1]; //
+ } //
+ }
+ else{
+ Msg::Warning("Cannot put a new point on Surface %d",gf->tag());
+ }
+ }
} //
} /// end non linear -------------------------------------------------//
@@ -325,6 +351,8 @@ bool compute4neighbors (GFace *gf, // the surface
void packingOfParallelograms(GFace* gf, std::vector<MVertex*> &packed, std::vector<SMetric3> &metrics){
#if defined(HAVE_RTREE)
+ const bool goNonLinear = true;
+
// FILE *f = fopen ("parallelograms.pos","w");
// get all the boundary vertices
@@ -347,7 +375,7 @@ void packingOfParallelograms(GFace* gf, std::vector<MVertex*> &packed, std::vec
SPoint2 newp[4][NUMDIR];
std::set<MVertex*>::iterator it = bnd_vertices.begin() ;
for (; it != bnd_vertices.end() ; ++it){
- compute4neighbors (gf, *it, false, newp, metricField);
+ compute4neighbors (gf, *it, goNonLinear, newp, metricField);
surfacePointWithExclusionRegion *sp =
new surfacePointWithExclusionRegion (*it, newp, metricField);
// fifo.push(sp);
@@ -381,7 +409,7 @@ void packingOfParallelograms(GFace* gf, std::vector<MVertex*> &packed, std::vec
GPoint gp = gf->point(parent->_p[i][dir]);
MFaceVertex *v = new MFaceVertex(gp.x(),gp.y(),gp.z(),gf,gp.u(),gp.v());
// printf(" %g %g %g %g\n",parent._center.x(),parent._center.y(),gp.u(),gp.v());
- compute4neighbors (gf, v, false, newp, metricField);
+ compute4neighbors (gf, v, goNonLinear, newp, metricField);
surfacePointWithExclusionRegion *sp =
new surfacePointWithExclusionRegion (v, newp, metricField, parent);
// fifo.push(sp);
diff --git a/benchmarks/hex/twoHoles.geo b/benchmarks/hex/twoHoles.geo
index a30c0806b5..f7f1c59894 100644
--- a/benchmarks/hex/twoHoles.geo
+++ b/benchmarks/hex/twoHoles.geo
@@ -3,7 +3,7 @@
Mesh.Algorithm = 9; //8 = delquad or 9= 2D R-tree
Mesh.Algorithm3D = 9; // 3D R-tree
Mesh.Recombine3DAll = 1;
-Mesh.RecombineAll = 1;
+//Mesh.RecombineAll = 1;
//Mesh.RecombinationAlgorithm = 2;
Mesh.Smoothing = 0;
--
GitLab