From 3f50aab22d5c9533fe011a01ed1dcb9b08f151e7 Mon Sep 17 00:00:00 2001
From: Jean-Francois Remacle <jean-francois.remacle@uclouvain.be>
Date: Wed, 21 Jan 2009 08:27:19 +0000
Subject: [PATCH] *** empty log message ***

---
 Common/OctreeInternals.cpp |  4 ++--
 Geo/MVertex.cpp            |  7 +++++--
 Mesh/Field.cpp             | 41 +++++++++++++++++++-------------------
 Mesh/meshGFace.cpp         |  2 ++
 Post/Makefile              |  2 +-
 Post/OctreePost.cpp        | 24 ++++++++++++++++++++++
 Post/OctreePost.h          |  2 ++
 Post/shapeFunctions.h      | 10 ++++++++--
 8 files changed, 65 insertions(+), 27 deletions(-)

diff --git a/Common/OctreeInternals.cpp b/Common/OctreeInternals.cpp
index f5ea8f5b83..600143b0a1 100644
--- a/Common/OctreeInternals.cpp
+++ b/Common/OctreeInternals.cpp
@@ -292,10 +292,10 @@ void *searchElement(octantBucket *_buckets_head, double *_pt, globalInfo *_globa
       flag = xyzInElement(ptrToEle, _pt);
     if (flag == 1) return ptrToEle;
   }
-        
+     
   ptrBucket = findElementBucket(_buckets_head, _pt);
   if (ptrBucket == NULL) {
-    // printf("Error! the point is not in the domain.\n");
+    printf("Error! the point is not in the domain.\n");
     return NULL;
   }     
 
diff --git a/Geo/MVertex.cpp b/Geo/MVertex.cpp
index 1a14720c41..46f399ed82 100644
--- a/Geo/MVertex.cpp
+++ b/Geo/MVertex.cpp
@@ -320,8 +320,11 @@ bool reparamMeshVertexOnFace(const MVertex *v, const GFace *gf, SPoint2 &param)
 
   if(v->onWhat()->dim() == 0){
     GVertex *gv = (GVertex*)v->onWhat();
-    param = gv->reparamOnFace(gf, 1);
-
+    // hack for bug in periodic curves
+    if (gv->getNativeType() == GEntity::GmshModel && gf->geomType() == GEntity::Plane)
+      param = gf->parFromPoint(SPoint3(v->x(), v->y(), v->z()));
+    else
+      param = gv->reparamOnFace(gf, 1);
     // shout, we could be on a seam
     std::list<GEdge*> ed = gv->edges();
     for(std::list<GEdge*>::iterator it = ed.begin(); it != ed.end(); it++)
diff --git a/Mesh/Field.cpp b/Mesh/Field.cpp
index 5267d91184..39bc5d4967 100644
--- a/Mesh/Field.cpp
+++ b/Mesh/Field.cpp
@@ -968,26 +968,27 @@ class PostViewField : public Field
     if(!octree->searchScalar(x, y, z, &l, 0)) {
       // try really hard to find an element around the point
       /*
-	double fact[4] = {1.e-6, 1.e-5, 1.e-4, 1.e-3};
-	for(int i = 0; i < 4; i++){
-	  double eps = CTX.lc * fact[i];
-	  printf("approx search witg eps=%g\n", eps);
-	  if(octree->searchScalar(x + eps, y, z, &l, 0)) break;
-	  if(octree->searchScalar(x - eps, y, z, &l, 0)) break;
-	  if(octree->searchScalar(x, y + eps, z, &l, 0)) break;
-	  if(octree->searchScalar(x, y - eps, z, &l, 0)) break;
-	  if(octree->searchScalar(x, y, z + eps, &l, 0)) break;
-	  if(octree->searchScalar(x, y, z - eps, &l, 0)) break;
-	  if(octree->searchScalar(x + eps, y - eps, z - eps, &l, 0)) break;
-	  if(octree->searchScalar(x + eps, y + eps, z - eps, &l, 0)) break;
-	  if(octree->searchScalar(x - eps, y - eps, z - eps, &l, 0)) break;
-	  if(octree->searchScalar(x - eps, y + eps, z - eps, &l, 0)) break;
-	  if(octree->searchScalar(x + eps, y - eps, z + eps, &l, 0)) break;
-	  if(octree->searchScalar(x + eps, y + eps, z + eps, &l, 0)) break;
-	  if(octree->searchScalar(x - eps, y - eps, z + eps, &l, 0)) break;
-	  if(octree->searchScalar(x - eps, y + eps, z + eps, &l, 0)) break;
-	}
-      */
+      double fact[4] = {1.e-6, 1.e-5, 1.e-4, 1.e-2};
+      for(int i = 0; i < 4; i++){
+	double eps = CTX.lc * fact[i];
+	//	printf("approx search witg eps=%g\n", eps);
+	if(octree->searchScalar(x + eps, y, z, &l, 0)) break;
+	if(octree->searchScalar(x - eps, y, z, &l, 0)) break;
+	if(octree->searchScalar(x, y + eps, z, &l, 0)) break;
+	if(octree->searchScalar(x, y - eps, z, &l, 0)) break;
+	if(octree->searchScalar(x, y, z + eps, &l, 0)) break;
+	if(octree->searchScalar(x, y, z - eps, &l, 0)) break;
+	if(octree->searchScalar(x + eps, y - eps, z - eps, &l, 0)) break;
+	if(octree->searchScalar(x + eps, y + eps, z - eps, &l, 0)) break;
+	if(octree->searchScalar(x - eps, y - eps, z - eps, &l, 0)) break;
+	if(octree->searchScalar(x - eps, y + eps, z - eps, &l, 0)) break;
+	if(octree->searchScalar(x + eps, y - eps, z + eps, &l, 0)) break;
+	if(octree->searchScalar(x + eps, y + eps, z + eps, &l, 0)) break;
+	if(octree->searchScalar(x - eps, y - eps, z + eps, &l, 0)) break;
+	if(octree->searchScalar(x - eps, y + eps, z + eps, &l, 0)) break;
+      } 
+      */     
+      //      printf("oops\n");
     }
     if(l <= 0 && crop_negative_values) return MAX_LC;
     return l;
diff --git a/Mesh/meshGFace.cpp b/Mesh/meshGFace.cpp
index 428f65fcfd..9657edecc1 100644
--- a/Mesh/meshGFace.cpp
+++ b/Mesh/meshGFace.cpp
@@ -432,6 +432,7 @@ static bool gmsh2DMeshGenerator(GFace *gf, int RECUR_ITER, bool debug = true)
     reparamMeshVertexOnFace(here, gf, param);
     U_[count] = param[0];
     V_[count] = param[1];
+    //    printf("coucou : %g %g -> %g %g\n",here->x(),here->y(),param.x(),param.y());
     (*itv)->setIndex(count);
     numbered_vertices[(*itv)->getIndex()] = *itv;
     bbox += SPoint3(param.x(), param.y(), 0);
@@ -1368,6 +1369,7 @@ void meshGFace::operator() (GFace *gf)
 
   Msg::Debug("Generating the mesh");
   if(noseam(gf) || gf->getNativeType() == GEntity::GmshModel || gf->edgeLoops.empty()){
+    //gmsh2DMeshGeneratorPeriodic(gf, debugSurface >= 0 || debugSurface == -100);
     gmsh2DMeshGenerator(gf, 0, debugSurface >= 0 || debugSurface == -100);
   }
   else{
diff --git a/Post/Makefile b/Post/Makefile
index 8e603f1eb5..1c5e550467 100644
--- a/Post/Makefile
+++ b/Post/Makefile
@@ -19,7 +19,7 @@ SRC = PView.cpp PViewIO.cpp\
           PViewDataList.cpp PViewDataListIO.cpp\
           PViewDataGModel.cpp PViewDataGModelIO.cpp\
         PViewOptions.cpp\
-      adaptiveData.cpp\
+      adaptiveData.cpp shapeFunctions.cpp\
       OctreePost.cpp\
       ColorTable.cpp
 
diff --git a/Post/OctreePost.cpp b/Post/OctreePost.cpp
index 3f70468da4..64b7c7d88c 100644
--- a/Post/OctreePost.cpp
+++ b/Post/OctreePost.cpp
@@ -20,6 +20,7 @@
 static void minmax(int n, double *X, double *Y, double *Z,
                    double *min, double *max)
 {
+  const double eps = 1.e-1;
   min[0] = X[0];
   min[1] = Y[0];
   min[2] = Z[0];
@@ -35,6 +36,15 @@ static void minmax(int n, double *X, double *Y, double *Z,
     max[1] = (Y[i] > max[1]) ? Y[i] : max[1];
     max[2] = (Z[i] > max[2]) ? Z[i] : max[2];
   }
+  const double L = eps * sqrt ((min[0]-max[0])*(min[0]-max[0])+
+			       (min[1]-max[1])*(min[1]-max[1])+
+			       (min[2]-max[2])*(min[2]-max[2]));
+  min[0] -= L;
+  min[1] -= L;
+  min[2] -= L;
+  max[0] += L;
+  max[1] += L;
+  max[2] += L;
 }
 
 static void centroid(int n, double *X, double *Y, double *Z, double *c)
@@ -385,6 +395,20 @@ bool OctreePost::_getValue(void *in, int nbComp, double P[3], int timestep,
 
 bool OctreePost::searchScalar(double x, double y, double z, double *values, 
                               int step, double *size)
+{
+  bool a =  searchScalar_(x,y,z,values,step,size);
+  if (!a){
+    element::setTolerance(10.);
+    a =  searchScalar_(x,y,z,values,step,size);
+    element::setTolerance(1.e-3);
+  }    
+  if (!a)printf("cannot find %g %g %g\n",x,y,z);
+  return a;
+}
+
+
+bool OctreePost::searchScalar_(double x, double y, double z, double *values, 
+                              int step, double *size)
 {
   double P[3] = {x, y, z};
 
diff --git a/Post/OctreePost.h b/Post/OctreePost.h
index dcf654c12a..172125d449 100644
--- a/Post/OctreePost.h
+++ b/Post/OctreePost.h
@@ -30,6 +30,8 @@ class OctreePost
                  double *elementSize);
   bool _getValue(void *in, int nbComp, double P[3], int step, 
                  double *values, double *elementSize);
+  bool searchScalar_(double x, double y, double z, double *values, 
+                    int step = -1, double *size = 0);
  public :
   OctreePost(PView *);
   ~OctreePost();
diff --git a/Post/shapeFunctions.h b/Post/shapeFunctions.h
index 5bb00835f8..da06e97770 100644
--- a/Post/shapeFunctions.h
+++ b/Post/shapeFunctions.h
@@ -9,15 +9,21 @@
 #include "Numeric.h"
 #include "GmshMessage.h"
 
-#define ONE     (1. + 1.e-6)
-#define ZERO    (-1.e-6)
 #define SQU(a)  ((a)*(a))
 
 class element{
 protected:
   double *_x, *_y, *_z;
+  static double ONE, ZERO;
 public:
   element(double *x, double *y, double *z) : _x(x), _y(y), _z(z) {}
+  static void setTolerance (const double tol){
+    ONE = 1.+tol;
+    ZERO = -tol;
+  }
+  static double getTolerance () {
+    return -ZERO;
+  }
   virtual ~element(){}
   virtual int getDimension() = 0;
   virtual int getNumNodes() = 0;
-- 
GitLab