Added some python access to adaptivity

Frontal mesher is less sensitive to small geometrical features (avoid obtuse triangles)

Geo/MVertex.h

@@ -170,4 +170,11 @@ bool reparamMeshVertexOnEdge(MVertex *v, const GEdge *ge, double &param);
 
 double angle3Vertices(MVertex *p1, MVertex *p2, MVertex *p3);
 
+inline double distance (MVertex *v1, MVertex *v2){
+  const double dx = v1->x() - v2->x();
+  const double dy = v1->y() - v2->y();
+  const double dz = v1->z() - v2->z();
+  return sqrt(dx*dx+dy*dy+dz*dz);
+}
+
 #endif

Mesh/meshGEdge.cpp

@@ -387,11 +387,12 @@ void meshGEdge::operator() (GEdge *ge)
   
   // force odd number of points for if blossom is used for recombination
   if(ge->meshAttributes.Method != MESH_TRANSFINITE &&
-     CTX::instance()->mesh.algoRecombine == 1 && N % 2 == 0){
-    if(CTX::instance()->mesh.recombineAll){
+     (CTX::instance()->mesh.algoRecombine == 1 || CTX::instance()->mesh.recombineAll) && N % 2 == 0){
+    //    if(CTX::instance()->mesh.recombineAll){
       N++;
+      //    }
   }
-    else{
+  /*    else{
       std::list<GFace*> faces = ge->faces();
       for(std::list<GFace*>::iterator it = faces.begin(); it != faces.end(); it++){
         if((*it)->meshAttributes.recombine){
@@ -401,6 +402,7 @@ void meshGEdge::operator() (GEdge *ge)
       }
     }
   }
+  */
 
   // printFandPrimitive(ge->tag(),Points);
 

Mesh/meshGFaceDelaunayInsertion.cpp

@@ -574,19 +574,16 @@ bool insertVertex(bool force, GFace *gf, MVertex *v, double *param , MTri3 *t,
     MTri3 *t4;
     t4 = new MTri3(t, LL,0,&Us,&Vs,gf);
 
-    double d1 = sqrt((it->v[0]->x() - v->x()) * (it->v[0]->x() - v->x()) +
-                     (it->v[0]->y() - v->y()) * (it->v[0]->y() - v->y()) +
-                     (it->v[0]->z() - v->z()) * (it->v[0]->z() - v->z()));
-    double d2 = sqrt((it->v[1]->x() - v->x()) * (it->v[1]->x() - v->x()) +
-                     (it->v[1]->y() - v->y()) * (it->v[1]->y() - v->y()) +
-                     (it->v[1]->z() - v->z()) * (it->v[1]->z() - v->z()));
-    const double MID[3] = {0.5*(it->v[0]->x()+it->v[1]->x()),
-                           0.5*(it->v[0]->y()+it->v[1]->y()),
-                           0.5*(it->v[0]->z()+it->v[1]->z())};
-    double d3 = sqrt((MID[0] - v->x()) * (MID[0] - v->x()) +
-                     (MID[1] - v->y()) * (MID[1] - v->y()) +
-                     (MID[2] - v->z()) * (MID[2] - v->z()));
-    if ((d1 < LL * .25 || d2 < LL * .25 || d3 < LL * .25) && !force) {
+    //    double din = t->getInnerRadius();
+
+    double d1 = distance(it->v[0],v);
+    double d2 = distance(it->v[1],v);
+    double d3 = distance(it->v[0],it->v[1]);
+
+    // avoid angles that are too obtuse
+    double cosv = ((d1*d1+d2*d2-d3*d3)/(2.*d1*d2));
+
+    if ((d1 < LL * .25 || d2 < LL * .25 || cosv < -0.5) && !force) {
       onePointIsTooClose = true;
     }
 
@@ -746,35 +743,8 @@ static bool insertAPoint(GFace *gf, std::set<MTri3*,compareTri3Ptr>::iterator it
   MTri3 *ptin = 0;
   bool inside = false;
   double uv[2];
-  // FIXME !!! ----> FIXED (JFR)
-  if (0 && MTri3::radiusNorm == 2){
-    inside = invMapUV(worst->tri(), center, Us, Vs, uv, 1.e-8);
-    if (inside)ptin = worst;
-    if (!inside && worst->getNeigh(0)){
-      inside |= invMapUV(worst->getNeigh(0)->tri(), center, Us, Vs, uv, 1.e-8);
-      if (inside)ptin = worst->getNeigh(0);
-    }
-    if (!inside && worst->getNeigh(1)){
-      inside |= invMapUV(worst->getNeigh(1)->tri(), center, Us, Vs, uv, 1.e-8);
-      if (inside)ptin = worst->getNeigh(1);
-    }
-    if (!inside && worst->getNeigh(2)){
-      inside |= invMapUV(worst->getNeigh(2)->tri(), center, Us, Vs, uv, 1.e-8);
-      if (inside)ptin = worst->getNeigh(2);
-    }
-    if (!inside){
-      ptin =  search4Triangle (worst, center, Us, Vs, AllTris,uv);
-      if (ptin)inside = true;
-    }
-  }
-  else {
   ptin = search4Triangle (worst, center, Us, Vs, AllTris,uv);
-    // if (!ptin){
-    //   printf("strange : %g %g seems to be out of the domain for face %d\n",
-    //          center[0],center[1],gf->tag());
-    // }
   if (ptin) inside = true;
-  }
 
   if (inside) {
     // we use here local coordinates as real coordinates
@@ -1025,24 +995,6 @@ double optimalPointFrontal (GFace *gf,
 			    2 * dir[1] * dir[0] * metric[1] +
 			    dir[1] * dir[1] * metric[2]);
 
-  //  const double p = 0.5 * lengthMetric(P, Q, metric); // / RATIO;
-  /*
-  const double rhoM1 = 0.5 *
-    (vSizes[base->getVertex(ip1)->getIndex()] +
-     vSizes[base->getVertex(ip2)->getIndex()] ) / sqrt(3.);// * RATIO;
-  const double rhoM2 = 0.5 *
-    (vSizesBGM[base->getVertex(ip1)->getIndex()] +
-     vSizesBGM[base->getVertex(ip2)->getIndex()] ) / sqrt(3.);// * RATIO;
-  const double rhoM  = Extend1dMeshIn2dSurfaces() ? std::min(rhoM1, rhoM2) : rhoM2;
-  */
-
-  //  const double rhoM_hat = std::min(std::max(rhoM, p), (p * p + q * q) / (2 * q));
-  //  double d = (rhoM_hat + sqrt (rhoM_hat * rhoM_hat - p * p)) / RATIO;
-
-  //const double d = 100./RATIO;
-  //  printf("(%g %g) (%g %g) %g %g %g %g %g %g\n",
-  //          P[0],P[1],Q[0],Q[1],RATIO,p,q,rhoM,rhoM_hat,d);
-  //  printf("size %12.5E\n",vSizesBGM[base->getVertex(ip1)->getIndex()]);
   const double rhoM1 = 0.5*
     (vSizes[base->getVertex(ip1)->getIndex()] +
      vSizes[base->getVertex(ip2)->getIndex()] ) ;// * RATIO;
@@ -1053,17 +1005,22 @@ double optimalPointFrontal (GFace *gf,
   const double rhoM_hat = rhoM;
 
   const double q = lengthMetric(center, midpoint, metric);
-
   const double d = rhoM_hat * sqrt(3.)*0.5;
 
+  // d is corrected in a way that the mesh size is computed at point newPoint
+  
+
   //  printf("%12.5E %12.5E\n",d,RATIO);
 
   //  const double L = d ;
   // avoid to go toooooo far
   const double L = d > q ? q : d;
   
+
   newPoint[0] = midpoint[0] + L * dir[0]/RATIO;
   newPoint[1] = midpoint[1] + L * dir[1]/RATIO;
+
+
   return L;// > q ? d : q;
 }
 

Mesh/meshGFaceOptimize.cpp

@@ -3039,6 +3039,8 @@ void recombineIntoQuads(GFace *gf,
     }
     double t2 = Cpu();
     Msg::Info("Blossom recombination algorithm completed (%g s)", t2 - t1);
+    quadsToTriangles(gf, 0.01);
+    //    if(haveParam) laplaceSmoothing(gf,CTX::instance()->mesh.nbSmoothing);
     return;
   }
 

Post/PViewAsSimpleFunction.cpp

@@ -12,6 +12,7 @@ double PViewEvaluator::operator() (const double x, const double y, const double
   PViewData * pvd = _pv->getData();
   double value;
   bool found = pvd->searchScalar(x, y, z, &value, _step);
+  printf("found %d %g %g %g %g\n",found,x,y,value,x*x+y*y);
   if (found) return value;
   return 1.e22;
 }

wrappers/gmshpy/gmshCommon.i

@@ -39,7 +39,6 @@ namespace std {
   %template(DoubleVector) std::vector<double, std::allocator<double> >;
   %template(DoubleVectorVector) std::vector<std::vector<double, std::allocator<double> > >;
   %template(StringVector) std::vector<std::string, std::allocator<std::string> >;
-  %template(IntDoubleMap) std::map < int , double >;
 }
 %pointer_functions(double,doublep)
 %pointer_functions(int,intp)

wrappers/gmshpy/gmshPost.i

@@ -2,6 +2,8 @@
 %module gmshPost
 %include typemaps.i
 %include std_string.i
+%include std_vector.i
+%include std_map.i
 
 %{
   #include "GmshConfig.h"
@@ -14,11 +16,20 @@
 #endif
 %}
 
+namespace std {
+  %template(DoubleVector) std::vector<double, std::allocator<double> >;
+  %template(IntDoubleVectorMap) std::map < int , std::vector<double, std::allocator<double> > >;
+}
+
+
 %include "GmshConfig.h"
 #if defined(HAVE_POST)
 %include "PView.h"
 %include "PViewFactory.h"
 %apply double &OUTPUT { double &val}
 %include "PViewData.h"
+%include "simpleFunction.h"
+%template(simpleFunctionDouble) simpleFunction<double>;
+%include "PViewAsSimpleFunction.h"
 #endif