Field.cpp

      if(!_f[iFunction])
        metr(index[iFunction][0], index[iFunction][1]) = MAX_LC;
      else{
	std::vector<double> values(3 + _fields[iFunction].size()), res(1);
	values[0] = x;
	values[1] = y;
	values[2] = z;
	int i = 3;
	for(std::set<int>::iterator it = _fields[iFunction].begin();
            it != _fields[iFunction].end(); it++){
	  Field *field = GModel::current()->getFields()->get(*it);
	  values[i++] = field ? (*field)(x, y, z) : MAX_LC;
	}
	if(_f[iFunction]->eval(values, res))
	  metr(index[iFunction][0],index[iFunction][1]) =  res[0];
	else
	  metr(index[iFunction][0],index[iFunction][1]) = MAX_LC;
      }
    }
  }
};

class MathEvalField : public Field
{
  MathEvalExpression expr;
  std::string f;

 public:
  void myAction(){
    printf("doing sthg \n");
  }
  MathEvalField()
  {
    options["F"] = new FieldOptionString
      (f, "Mathematical function to evaluate.", &update_needed);
    f = "F2 + Sin(z)";
    callbacks["test"] = new FieldCallbackGeneric<MathEvalField>(this, &MathEvalField::myAction, "description blabla");
  }
  double operator() (double x, double y, double z, GEntity *ge=0)
  {
    if(update_needed) {
      if(!expr.set_function(f))
        Msg::Error("Field %i: Invalid matheval expression \"%s\"",
                   this->id, f.c_str());
      update_needed = false;
    }
    return expr.evaluate(x, y, z);
  }
  const char *getName()
  {
    return "MathEval";
  }
  std::string getDescription()
  {
    return "Evaluate a mathematical expression. The expression can contain "
      "x, y, z for spatial coordinates, F0, F1, ... for field values, and "
      "and mathematical functions.";
  }
};

class MathEvalFieldAniso : public Field
{
  MathEvalExpressionAniso expr;
  std::string f[6];
 public:
  virtual bool isotropic () const { return false; }
  MathEvalFieldAniso()
  {
    options["m11"] = new FieldOptionString
      (f[0], "element 11 of the metric tensor.", &update_needed);
    f[0] = "F2 + Sin(z)";
    options["m22"] = new FieldOptionString
      (f[1], "element 22 of the metric tensor.", &update_needed);
    f[1] = "F2 + Sin(z)";
    options["m33"] = new FieldOptionString
      (f[2], "element 33 of the metric tensor.", &update_needed);
    f[2] = "F2 + Sin(z)";
    options["m12"] = new FieldOptionString
      (f[3], "element 12 of the metric tensor.", &update_needed);
    f[3] = "F2 + Sin(z)";
    options["m13"] = new FieldOptionString
      (f[4], "element 13 of the metric tensor.", &update_needed);
    f[4] = "F2 + Sin(z)";
    options["m23"] = new FieldOptionString
      (f[5], "element 23 of the metric tensor.", &update_needed);
    f[5] = "F2 + Sin(z)";
  }
  void operator() (double x, double y, double z, SMetric3 &metr, GEntity *ge=0)
  {
    if(update_needed) {
      for (int i=0;i<6;i++){
	if(!expr.set_function(i,f[i]))
	  Msg::Error("Field %i: Invalid matheval expression \"%s\"",
		     this->id, f[i].c_str());
      }
      update_needed = false;
    }
    expr.evaluate(x, y, z, metr);
  }
  double operator() (double x, double y, double z, GEntity *ge=0)
  {
    if(update_needed) {
      for (int i = 0; i < 6; i++){
	if(!expr.set_function(i, f[i]))
	  Msg::Error("Field %i: Invalid matheval expression \"%s\"",
		     this->id, f[i].c_str());
      }
      update_needed = false;
    }
    SMetric3 metr;
    expr.evaluate(x, y, z, metr);
    return metr(0, 0);
  }
  const char *getName()
  {
    return "MathEvalAniso";
  }
  std::string getDescription()
  {
    return "Evaluate a metric expression. The expressions can contain "
      "x, y, z for spatial coordinates, F0, F1, ... for field values, and "
      "and mathematical functions.";
  }
};

class ParametricField : public Field
{
  MathEvalExpression expr[3];
  std::string f[3];
  int iField;
 public:
  ParametricField()
  {
    iField = 1;
    options["IField"] = new FieldOptionInt
      (iField, "Field index");
    options["FX"] = new FieldOptionString
      (f[0], "X component of parametric function", &update_needed);
    options["FY"] = new FieldOptionString
      (f[1], "Y component of parametric function", &update_needed);
    options["FZ"] = new FieldOptionString
      (f[2], "Z component of parametric function", &update_needed);
  }
  std::string getDescription()
  {
    return "Evaluate Field IField in parametric coordinates:\n\n"
      "  F = Field[IField](FX,FY,FZ)\n\n"
      "See the MathEval Field help to get a description of valid FX, FY "
      "and FZ expressions.";
  }
  double operator() (double x, double y, double z, GEntity *ge=0)
  {
    if(update_needed) {
      for(int i = 0; i < 3; i++) {
        if(!expr[i].set_function(f[i]))
          Msg::Error("Field %i : Invalid matheval expression \"%s\"",
                     this->id, f[i].c_str());
      }
      update_needed = false;
    }
    Field *field = GModel::current()->getFields()->get(iField);
    if(!field || iField == id) return MAX_LC;
    return (*field)(expr[0].evaluate(x, y, z),
                    expr[1].evaluate(x, y, z),
                    expr[2].evaluate(x, y, z));
  }
  const char *getName()
  {
    return "Param";
  }
};

#if defined(HAVE_POST)

class PostViewField : public Field
{
  OctreePost *octree;
  int view_index;
  bool crop_negative_values;
 public:
  PostViewField()
  {
    octree = 0;
    view_index = 0;
    options["IView"] = new FieldOptionInt
      (view_index, "Post-processing view index", &update_needed);
    crop_negative_values = true;
    options["CropNegativeValues"] = new FieldOptionBool
      (crop_negative_values, "return LC_MAX instead of a negative value (this "
       "option is needed for backward compatibility with the BackgroundMesh option");
  }
  ~PostViewField()
  {
    if(octree) delete octree;
  }
  PView *getView() const
  {
    // we should maybe test the unique view num instead, but that
    // would be slower
    if(view_index < 0 || view_index >= (int)PView::list.size()){
      Msg::Error("View[%d] does not exist", view_index);
      return 0;
    }
    PView *v = PView::list[view_index];
    if(v->getData()->hasModel(GModel::current())){
      Msg::Error("Cannot use view based on current mesh for background mesh: you might"
                 " want to use a list-based view (.pos file) instead");
      return 0;
    }
    return v;
  }
  virtual bool isotropic () const
  {
    PView *v = getView();
    if(v && v->getData()->getNumTensors()) return false;
    return true;
  }
  double operator() (double x, double y, double z, GEntity *ge=0)
  {
    PView *v = getView();
    if(!v) return MAX_LC;
    if(update_needed){
      if(octree) delete octree;
      octree = new OctreePost(v);
      update_needed = false;
    }
    double l = 0.;
    // use large tolerance (in element reference coordinates) to maximize chance
    // of finding an element
    if(!octree->searchScalarWithTol(x, y, z, &l, 0, 0, 1.))
      Msg::Info("No scalar element found containing point (%g,%g,%g)", x, y, z);
    if(l <= 0 && crop_negative_values) return MAX_LC;
    return l;
  }
  void operator() (double x, double y, double z, SMetric3 &metr, GEntity *ge=0)
  {
    PView *v = getView();
    if(!v) return;
    if(update_needed){
      if(octree) delete octree;
      octree = new OctreePost(v);
      update_needed = false;
    }
    double l[9];
    // use large tolerance (in element reference coordinates) to maximize chance
    // of finding an element
    if(!octree->searchTensorWithTol(x, y, z, l, 0, 0, 1.)){
      Msg::Info("No tensor element found containing point (%g,%g,%g)", x, y, z);
      return;
    }
    metr(0, 0) = l[0];
    metr(0, 1) = l[1];
    metr(0, 2) = l[2];
    metr(1, 0) = l[3];
    metr(1, 1) = l[4];
    metr(1, 2) = l[5];
    metr(2, 0) = l[6];
    metr(2, 1) = l[7];
    metr(2, 2) = l[8];
  }
  const char *getName()
  {
    return "PostView";
  }
  std::string getDescription()
  {
    return "Evaluate the post processing view IView.";
  }
};

#endif

class MinAnisoField : public Field
{
  std::list<int> idlist;
 public:
  MinAnisoField()
  {
    options["FieldsList"] = new FieldOptionList
      (idlist, "Field indices", &update_needed);
  }
  virtual bool isotropic () const {return false;}
  std::string getDescription()
  {
    return "Take the intersection of a list of possibly anisotropic fields.";
  }
  virtual void operator() (double x, double y, double z, SMetric3 &metr, GEntity *ge=0)
  {
    SMetric3 v (1./MAX_LC);
    for(std::list<int>::iterator it = idlist.begin(); it != idlist.end(); it++) {
      Field *f = (GModel::current()->getFields()->get(*it));
      SMetric3 ff;
      if(f && *it != id) {
	if (f->isotropic()){
	  double l = (*f) (x, y, z, ge);
	  ff = SMetric3(1./(l*l));
	}
	else{
	  (*f) (x, y, z, ff, ge);
	}
	v = intersection_conserve_mostaniso(v,ff);
      }
    }
    metr = v;
  }
  double operator() (double x, double y, double z, GEntity *ge=0)
  {
    SMetric3 metr (1./MAX_LC);
    double v = MAX_LC;
    for(std::list<int>::iterator it = idlist.begin(); it != idlist.end(); it++) {
      Field *f = (GModel::current()->getFields()->get(*it));
      SMetric3 m;
      if(f && *it != id){
        if (!f->isotropic()){
          (*f)(x, y, z, m, ge);
        }
	else {
          double L = (*f)(x, y, z, ge);
          for (int i = 0; i < 3; i++)
            m(i,i) = 1. / (L*L);
	}
      }
      metr = intersection(metr, m);
    }
    fullMatrix<double> V(3,3);
    fullVector<double> S(3);
    metr.eig(V, S, 1);
    double val = sqrt(1./S(2)); //S(2) is largest eigenvalue
    return std::min(v, val);
  }
  const char *getName()
  {
    return "MinAniso";
  }
};

class MinField : public Field
{
  std::list<int> idlist;
 public:
  MinField()
  {
    options["FieldsList"] = new FieldOptionList
      (idlist, "Field indices", &update_needed);
  }
  std::string getDescription()
  {
    return "Take the minimum value of a list of fields.";
  }
  double operator() (double x, double y, double z, GEntity *ge=0)
  {
    double v = MAX_LC;
    for(std::list<int>::iterator it = idlist.begin(); it != idlist.end(); it++) {
      Field *f = (GModel::current()->getFields()->get(*it));
      if(f && *it != id) v = std::min(v, (*f) (x, y, z, ge));
    }
    return v;
  }
  const char *getName()
  {
    return "Min";
  }
};

class MaxField : public Field
{
  std::list<int> idlist;
 public:
  MaxField()
  {
    options["FieldsList"] = new FieldOptionList
      (idlist, "Field indices", &update_needed);
  }
  std::string getDescription()
  {
    return "Take the maximum value of a list of fields.";
  }
  double operator() (double x, double y, double z, GEntity *ge=0)
  {
    double v = -MAX_LC;
    for(std::list<int>::iterator it = idlist.begin(); it != idlist.end(); it++) {
      Field *f = (GModel::current()->getFields()->get(*it));
      if(f && *it != id) v = std::max(v, (*f) (x, y, z, ge));
    }
    return v;
  }
  const char *getName()
  {
    return "Max";
  }
};

class RestrictField : public Field
{
  int iField;
  std::list<int> edges, faces, regions;
 public:
  RestrictField()
  {
    iField = 1;
    options["IField"] = new FieldOptionInt(iField, "Field index");
    options["EdgesList"] = new FieldOptionList(edges, "Curve indices");
    options["FacesList"] = new FieldOptionList(faces, "Surface indices");
    options["RegionsList"] = new FieldOptionList(regions, "Volume indices");
  }
  std::string getDescription()
  {
    return "Restrict the application of a field to a given list of geometrical "
      "curves, surfaces or volumes.";
  }
  double operator() (double x, double y, double z, GEntity *ge=0)
  {
    Field *f = (GModel::current()->getFields()->get(iField));
    if(!f || iField == id) return MAX_LC;
    if(!ge) return (*f) (x, y, z);
    if((ge->dim() == 0) ||
       (ge->dim() == 1 && std::find
        (edges.begin(), edges.end(), ge->tag()) != edges.end()) ||
       (ge->dim() == 2 && std::find
        (faces.begin(), faces.end(), ge->tag()) != faces.end()) ||
       (ge->dim() == 3 && std::find
        (regions.begin(), regions.end(), ge->tag()) != regions.end()))
      return (*f) (x, y, z);
    return MAX_LC;
  }
  const char *getName()
  {
    return "Restrict";
  }
};

#if defined(HAVE_ANN)
struct AttractorInfo{
  AttractorInfo (int a=0, int b=0, double c=0, double d=0)
    : ent(a),dim(b),u(c),v(d) {
  }
  int ent,dim;
  double u,v;
};

class AttractorAnisoCurveField : public Field {
  ANNkd_tree *kdtree;
  ANNpointArray zeronodes;
  ANNidxArray index;
  ANNdistArray dist;
  std::list<int> edges_id;
  double dMin, dMax, lMinTangent, lMaxTangent, lMinNormal, lMaxNormal;
  int n_nodes_by_edge;
  std::vector<SVector3> tg;
  public:
  AttractorAnisoCurveField() : kdtree(0), zeronodes(0)
  {
    index = new ANNidx[1];
    dist = new ANNdist[1];
    n_nodes_by_edge = 20;
    update_needed = true;
    dMin = 0.1;
    dMax = 0.5;
    lMinNormal = 0.05;
    lMinTangent = 0.5;
    lMaxNormal = 0.5;
    lMaxTangent = 0.5;
    options["EdgesList"] = new FieldOptionList
      (edges_id, "Indices of curves in the geometric model", &update_needed);
    options["NNodesByEdge"] = new FieldOptionInt
      (n_nodes_by_edge, "Number of nodes used to discretized each curve",
       &update_needed);
    options["dMin"] = new FieldOptionDouble
      (dMin, "Minimum distance, bellow this distance from the curves, "
       "prescribe the minimum mesh sizes.");
    options["dMax"] = new FieldOptionDouble
      (dMax, "Maxmium distance, above this distance from the curves, prescribe "
       "the maximum mesh sizes.");
    options["lMinTangent"] = new FieldOptionDouble
      (lMinTangent, "Minimum mesh size in the direction tangeant to the closest curve.");
    options["lMaxTangent"] = new FieldOptionDouble
      (lMaxTangent, "Maximum mesh size in the direction tangeant to the closest curve.");
    options["lMinNormal"] = new FieldOptionDouble
      (lMinNormal, "Minimum mesh size in the direction normal to the closest curve.");
    options["lMaxNormal"] = new FieldOptionDouble
      (lMaxNormal, "Maximum mesh size in the direction normal to the closest curve.");
  }
  virtual bool isotropic () const {return false;}
  ~AttractorAnisoCurveField()
  {
    if(kdtree) delete kdtree;
    if(zeronodes) annDeallocPts(zeronodes);
    delete[]index;
    delete[]dist;
  }
  const char *getName()
  {
    return "AttractorAnisoCurve";
  }
  std::string getDescription()
  {
    return "Compute the distance from the nearest curve in a list. Then the mesh "
           "size can be specified independently in the direction normal to the curve "
           "and in the direction parallel to the curve (Each curve "
           "is replaced by NNodesByEdge equidistant nodes and the distance from those "
           "nodes is computed.)";
  }
  void update()
  {
    if(zeronodes) {
      annDeallocPts(zeronodes);
      delete kdtree;
    }
    int totpoints = n_nodes_by_edge * edges_id.size();
    if(totpoints){
      zeronodes = annAllocPts(totpoints, 3);
    }
    tg.resize(totpoints);
    int k = 0;
    for(std::list<int>::iterator it = edges_id.begin();
        it != edges_id.end(); ++it) {
      Curve *c = FindCurve(*it);
      if(c) {
        for(int i = 1; i < n_nodes_by_edge-1; i++) {
          double u = (double)i / (n_nodes_by_edge - 1);
          Vertex V = InterpolateCurve(c, u, 0);
          zeronodes[k][0] = V.Pos.X;
          zeronodes[k][1] = V.Pos.Y;
          zeronodes[k][2] = V.Pos.Z;
          Vertex V2 = InterpolateCurve(c, u, 1);
          tg[k] = SVector3(V2.Pos.X, V2.Pos.Y, V2.Pos.Z);
          tg[k].normalize();
          k++;
        }
      }
      else {
        GEdge *e = GModel::current()->getEdgeByTag(*it);
        if(e) {
          for(int i = 1; i < n_nodes_by_edge-1; i++) {
            double u = (double)i / (n_nodes_by_edge - 1);
            Range<double> b = e->parBounds(0);
            double t = b.low() + u * (b.high() - b.low());
            GPoint gp = e->point(t);
            SVector3 d = e->firstDer(t);
            zeronodes[k][0] = gp.x();
            zeronodes[k][1] = gp.y();
            zeronodes[k][2] = gp.z();
            tg[k] = d;
            tg[k].normalize();
            k++;
          }
        }
      }
    }
    kdtree = new ANNkd_tree(zeronodes, totpoints, 3);
    update_needed = false;
  }
  void operator() (double x, double y, double z, SMetric3 &metr, GEntity *ge=0)
  {
    if(update_needed)
      update();
    double xyz[3] = { x, y, z };
    kdtree->annkSearch(xyz, 1, index, dist);
    double d = sqrt(dist[0]);
    double lTg = d < dMin ? lMinTangent : d > dMax ? lMaxTangent :
      lMinTangent + (lMaxTangent - lMinTangent) * (d - dMin) / (dMax - dMin);
    double lN = d < dMin ? lMinNormal : d > dMax ? lMaxNormal :
      lMinNormal + (lMaxNormal - lMinNormal) * (d - dMin) / (dMax - dMin);
    SVector3 t = tg[index[0]];
    SVector3 n0 = crossprod(t, fabs(t(0)) > fabs(t(1)) ? SVector3(0,1,0) :
                            SVector3(1,0,0));
    SVector3 n1 = crossprod(t, n0);
    metr = SMetric3(1/lTg/lTg, 1/lN/lN, 1/lN/lN, t, n0, n1);
  }
  virtual double operator() (double X, double Y, double Z, GEntity *ge=0)
  {
    if(update_needed)
      update();
    double xyz[3] = { X, Y, Z };
    kdtree->annkSearch(xyz, 1, index, dist);
    double d = sqrt(dist[0]);
    return std::max(d, 0.05);
  }
};

class AttractorField : public Field
{
  ANNkd_tree *kdtree;
  ANNpointArray zeronodes;
  ANNidxArray index;
  ANNdistArray dist;
  std::list<int> nodes_id, edges_id, faces_id;
  std::vector<AttractorInfo> _infos;
  int _xFieldId, _yFieldId, _zFieldId;
  Field *_xField, *_yField, *_zField;
  int n_nodes_by_edge;
 public:
  AttractorField(int dim, int tag, int nbe)
    : kdtree(0), zeronodes(0), n_nodes_by_edge(nbe)
  {
    index = new ANNidx[1];
    dist = new ANNdist[1];
    if (dim == 0) nodes_id.push_back(tag);
    else if (dim == 1) edges_id.push_back(tag);
    else if (dim == 2) faces_id.push_back(tag);
    _xFieldId = _yFieldId = _zFieldId = -1;
    update_needed = true;
  }
  AttractorField() : kdtree(0), zeronodes(0)
  {
    index = new ANNidx[1];
    dist = new ANNdist[1];
    n_nodes_by_edge = 20;
    options["NodesList"] = new FieldOptionList
      (nodes_id, "Indices of nodes in the geometric model", &update_needed);
    options["EdgesList"] = new FieldOptionList
      (edges_id, "Indices of curves in the geometric model", &update_needed);
    options["NNodesByEdge"] = new FieldOptionInt
      (n_nodes_by_edge, "Number of nodes used to discretized each curve",
       &update_needed);
    options["FacesList"] = new FieldOptionList
      (faces_id, "Indices of surfaces in the geometric model (Warning, this feature "
       "is still experimental. It might (read: will probably) give wrong results "
       "for complex surfaces)", &update_needed);
    _xFieldId = _yFieldId = _zFieldId = -1;
    options["FieldX"] = new FieldOptionInt
      (_xFieldId, "Id of the field to use as x coordinate.", &update_needed);
    options["FieldY"] = new FieldOptionInt
      (_yFieldId, "Id of the field to use as y coordinate.", &update_needed);
    options["FieldZ"] = new FieldOptionInt
      (_zFieldId, "Id of the field to use as z coordinate.", &update_needed);
  }
  ~AttractorField()
  {
    if(kdtree) delete kdtree;
    if(zeronodes) annDeallocPts(zeronodes);
    delete[]index;
    delete[]dist;
  }
  const char *getName()
  {
    return "Attractor";
  }
  std::string getDescription()
  {
    return "Compute the distance from the nearest node in a list. It can also "
      "be used to compute the distance from curves, in which case each curve "
      "is replaced by NNodesByEdge equidistant nodes and the distance from those "
      "nodes is computed.";
  }
  void getCoord(double x, double y, double z, double &cx, double &cy, double &cz,
                GEntity *ge = NULL) {
    cx = _xField ? (*_xField)(x, y, z, ge) : x;
    cy = _yField ? (*_yField)(x, y, z, ge) : y;
    cz = _zField ? (*_zField)(x, y, z, ge) : z;
  }
  std::pair<AttractorInfo,SPoint3> getAttractorInfo() const
  {
    return std::make_pair(_infos[index[0]], SPoint3(zeronodes[index[0]][0],
                                                    zeronodes[index[0]][1],
                                                    zeronodes[index[0]][2]));
  }
  virtual double operator() (double X, double Y, double Z, GEntity *ge=0)
  {
    _xField = _xFieldId >= 0 ? (GModel::current()->getFields()->get(_xFieldId)) : NULL;
    _yField = _yFieldId >= 0 ? (GModel::current()->getFields()->get(_yFieldId)) : NULL;
    _zField = _zFieldId >= 0 ? (GModel::current()->getFields()->get(_zFieldId)) : NULL;

    if(update_needed) {
      if(zeronodes) {
        annDeallocPts(zeronodes);
        delete kdtree;
      }

      std::vector<SPoint3> points;
      std::vector<SPoint2> uvpoints;
      std::vector<int> offset;
      offset.push_back(0);
      for(std::list<int>::iterator it = faces_id.begin();
          it != faces_id.end(); ++it) {
	GFace *f = GModel::current()->getFaceByTag(*it);
	if (f){

	  if (f->mesh_vertices.size()){
	    for (unsigned int i=0;i<f->mesh_vertices.size();i++){
	      MVertex *v = f->mesh_vertices[i];
	      double uu,vv;
	      v->getParameter(0,uu);
	      v->getParameter(1,vv);
	      points.push_back(SPoint3(v->x(),v->y(),v->z()));
	      uvpoints.push_back(SPoint2(uu,vv));
	    }
	  }
	  else {
	    SBoundingBox3d bb = f->bounds();
	    SVector3 dd = bb.max() - bb.min();
	    double maxDist = dd.norm() / n_nodes_by_edge ;
	    f->fillPointCloud(maxDist, &points, &uvpoints);
	    offset.push_back(points.size());
	  }
	}
      }

      int totpoints =
	nodes_id.size() +
	(n_nodes_by_edge-2) * edges_id.size() +
        ((points.size()) ? points.size() :
         n_nodes_by_edge * n_nodes_by_edge * faces_id.size());

      Msg::Info("%d points found in points clouds (%d edges)", totpoints,
                (int)edges_id.size());

      if(totpoints){
        zeronodes = annAllocPts(totpoints, 3);
        _infos.resize(totpoints);
      }
      int k = 0;
      for(std::list<int>::iterator it = nodes_id.begin();
          it != nodes_id.end(); ++it) {
	GVertex *gv = GModel::current()->getVertexByTag(*it);
	if(gv) {
	  getCoord(gv->x(), gv->y(), gv->z(), zeronodes[k][0],
		   zeronodes[k][1], zeronodes[k][2], gv);
	  _infos[k++] = AttractorInfo(*it,0,0,0);
        }
      }
      for(std::list<int>::iterator it = edges_id.begin();
          it != edges_id.end(); ++it) {
	GEdge *e = GModel::current()->getEdgeByTag(*it);
	if(e) {
	  if (e->mesh_vertices.size()){
	    for(unsigned int i = 0; i < e->mesh_vertices.size(); i++) {
	      double u ; e->mesh_vertices[i]->getParameter(0,u);
	      GPoint gp = e->point(u);
	      getCoord(gp.x(), gp.y(), gp.z(), zeronodes[k][0],
		       zeronodes[k][1], zeronodes[k][2], e);
	      _infos[k++] = AttractorInfo(*it,1,u,0);
	    }
	  }
	  int NNN = n_nodes_by_edge - e->mesh_vertices.size();
	  for(int i = 1; i < NNN - 1; i++) {
	    double u = (double)i / (NNN - 1);
	    Range<double> b = e->parBounds(0);
	    double t = b.low() + u * (b.high() - b.low());
	    GPoint gp = e->point(t);
	    getCoord(gp.x(), gp.y(), gp.z(), zeronodes[k][0],
		     zeronodes[k][1], zeronodes[k][2], e);
	    _infos[k++] = AttractorInfo(*it,1,t,0);
          }
        }
      }
      // This can lead to weird results as we generate attractors over
      // the whole parametric plane (we should really use a mesh,
      // e.g. a refined STL.)
      int count = 0;
      for(std::list<int>::iterator it = faces_id.begin();
          it != faces_id.end(); ++it) {
	GFace *f = GModel::current()->getFaceByTag(*it);
	if(f) {
	  if (points.size()){
	    for(int j = offset[count]; j < offset[count+1];j++) {
	      zeronodes[k][0] = points[j].x();
	      zeronodes[k][1] = points[j].y();
	      zeronodes[k][2] = points[j].z();
	      _infos[k++] = AttractorInfo(*it,2,uvpoints[j].x(),uvpoints[j].y());
	    }
	    count++;
	  }
	  else{
	    for(int i = 0; i < n_nodes_by_edge; i++) {
	      for(int j = 0; j < n_nodes_by_edge; j++) {
		double u = (double)i / (n_nodes_by_edge - 1);
		double v = (double)j / (n_nodes_by_edge - 1);
		Range<double> b1 = f->parBounds(0);
		Range<double> b2 = f->parBounds(1);
		double t1 = b1.low() + u * (b1.high() - b1.low());
		double t2 = b2.low() + v * (b2.high() - b2.low());
		GPoint gp = f->point(t1, t2);
		getCoord(gp.x(), gp.y(), gp.z(), zeronodes[k][0],
			 zeronodes[k][1], zeronodes[k][2], f);
		_infos[k++] = AttractorInfo(*it,2,u,v);
	      }
	    }
	  }
	}
	else {
	  printf("face %d not yet created\n",*it);
	}
      }
      kdtree = new ANNkd_tree(zeronodes, totpoints, 3);
      update_needed = false;
    }
    double xyz[3];
    getCoord(X, Y, Z, xyz[0], xyz[1], xyz[2], ge);
    kdtree->annkSearch(xyz, 1, index, dist);
    return sqrt(dist[0]);
  }
};

const char *BoundaryLayerField::getName()
{
  return "BoundaryLayer";
}

std::string BoundaryLayerField::getDescription()
{
  return "hwall * ratio^(dist/hwall)";
}

BoundaryLayerField::BoundaryLayerField()
{
  hwall_n = .1;
  hwall_t = .5;
  hfar = 1;
  ratio = 1.1;
  thickness = 1.e-2;
  fan_angle = 30;
  tgt_aniso_ratio = 1.e10;
  iRecombine = 0;
  iIntersect = 0;
  options["NodesList"] = new FieldOptionList
    (nodes_id, "Indices of nodes in the geometric model", &update_needed);
  options["EdgesList"] = new FieldOptionList
    (edges_id, "Indices of curves in the geometric model for which a boundary "
     "layer is needed", &update_needed);
  options["FacesList"] = new FieldOptionList
    (faces_id, "Indices of faces in the geometric model for which a boundary "
     "layer is needed", &update_needed);
  options["Quads"] = new FieldOptionInt
    (iRecombine, "Generate recombined elements in the boundary layer");
  options["IntersectMetrics"] = new FieldOptionInt
    (iIntersect, "Intersect metrics of all faces");
  options["hwall_n"] = new FieldOptionDouble
    (hwall_n, "Mesh Size Normal to the The Wall");
  options["fan_angle"] = new FieldOptionDouble
    (fan_angle, "Threshold angle for creating a mesh fan in the boundary layer");
  options["AnisoMax"] = new FieldOptionDouble
    (tgt_aniso_ratio, "Threshold angle for creating a mesh fan in the boundary layer");
  options["hwall_t"] = new FieldOptionDouble
    (hwall_t, "Mesh Size Tangent to the Wall");
  options["ratio"] = new FieldOptionDouble
    (ratio, "Size Ratio Between Two Successive Layers");
  options["hfar"] = new FieldOptionDouble
    (hfar, "Element size far from the wall");
  options["thickness"] = new FieldOptionDouble
    (thickness, "Maximal thickness of the boundary layer");
}

void BoundaryLayerField::removeAttractors()
{
  for (std::list<AttractorField *>::iterator it =  _att_fields.begin();
       it !=  _att_fields.end() ; ++it) delete *it;
  _att_fields.clear();
  update_needed = true;
}

void BoundaryLayerField::setupFor2d(int iF)
{
  if (1 || faces_id.size()){
    /* preprocess data in the following way
       remove GFaces from the attarctors (only used in 2D)
       for edges and vertices
     */
    if ( !faces_id_saved.size()){
      faces_id_saved = faces_id;
      edges_id_saved = edges_id;
      nodes_id_saved = nodes_id;
    }
    nodes_id.clear();
    edges_id.clear();
    faces_id.clear();

    //    printf("have %d %d\n",faces_id_saved.size(),edges_id_saved.size());

    ///  FIXME :
    /// NOT REALLY A NICE WAY TO DO IT (VERY AD HOC)
    /// THIS COULD BE PART OF THE INPUT
    /// OR (better) CHANGE THE PHILOSOPHY

    GFace *gf = GModel::current()->getFaceByTag(iF);
    std::list<GEdge*> ed = gf->edges();
    //    printf("face %d has %d edges\n",iF,ed.size());
    for (std::list<GEdge*>::iterator it = ed.begin();
	 it != ed.end() ; ++it){
      bool isIn = false;
      int iE = (*it)->tag();
      bool found = std::find ( edges_id_saved.begin(),edges_id_saved.end(),iE ) != edges_id_saved.end();
      //      printf("edges %d found %d\n",iE,found);
      // this edge is a BL Edge
      if (found){
	std::list<GFace*> fc = (*it)->faces();
	// one only face --> 2D --> BL
	if (fc.size() == 1) isIn = true;
	else {
	  // more than one face and
	  std::list<GFace*>::iterator itf = fc.begin();
	  bool found_this = std::find ( faces_id_saved.begin(),faces_id_saved.end(),iF ) != faces_id_saved.end();
	  if (!found_this)isIn = true;
	  else {
	    bool foundAll = true;
	    for ( ; itf != fc.end() ; ++itf){
	      int iF2 = (*itf)->tag();
	      foundAll &= std::find ( faces_id_saved.begin(),faces_id_saved.end(),iF2 ) != faces_id_saved.end();
	    }
	    if (foundAll)isIn = true;
	  }
	}
      }
      if (isIn){
	edges_id.push_back(iE);
	nodes_id.push_back ((*it)->getBeginVertex()->tag());
	nodes_id.push_back ((*it)->getEndVertex()->tag());
      }
    }
    printf("face %d %d BL Edges\n", iF, (int)edges_id.size());

    removeAttractors();
  }
}

void BoundaryLayerField::setupFor3d()
{
  faces_id = faces_id_saved;
  removeAttractors();
}


double BoundaryLayerField::operator() (double x, double y, double z, GEntity *ge)
{

  if (update_needed){
    for(std::list<int>::iterator it = nodes_id.begin();
	it != nodes_id.end(); ++it) {
      _att_fields.push_back(new AttractorField(0,*it,100000));
    }
    for(std::list<int>::iterator it = edges_id.begin();
	it != edges_id.end(); ++it) {
      _att_fields.push_back(new AttractorField(1,*it,300000));
    }
    for(std::list<int>::iterator it = faces_id.begin();
	it != faces_id.end(); ++it) {
      _att_fields.push_back(new AttractorField(2,*it,1200));
    }
    update_needed = false;
  }

  double dist = 1.e22;
  //AttractorField *cc;
  for (std::list<AttractorField*>::iterator it = _att_fields.begin();
       it != _att_fields.end(); ++it){
    double cdist = (*(*it)) (x, y, z);
    if (cdist < dist){
      //cc = *it;
      dist = cdist;
    }
  }
  current_distance = dist;
  //  const double dist = (*field) (x, y, z);
  //  current_distance = dist;
  const double lc = dist*(ratio-1) + hwall_t;
  //    double lc = hwall * pow (ratio, dist / hwall);
  return std::min (hfar,lc);
}

// assume that the closest point is one of the model vertices
void BoundaryLayerField::computeFor1dMesh (double x, double y, double z,
					   SMetric3 &metr)
{
  double xpk = 0., ypk = 0., zpk = 0.;
  double distk = 1.e22;
  for(std::list<int>::iterator it = nodes_id.begin();
      it != nodes_id.end(); ++it) {
    GVertex *v = GModel::current()->getVertexByTag(*it);
    double xp = v->x();
    double yp = v->y();
    double zp = v->z();
    const double dist = sqrt ((x - xp) *(x - xp)+
			      (y - yp) *(y - yp)+
			      (z - zp) *(z - zp));
    if (dist < distk){
      distk = dist;
      xpk = xp;
      ypk = yp;
      zpk = zp;
    }
  }

  const double ll1   = (distk*(ratio-1) + hwall_n) / (1. + 0.5 * (ratio - 1));
  // const double ll1   = (distk*(ratio-1) + hwall_n) / (1.);
  double lc_n  = std::min(ll1,hfar);

  if (distk > thickness) lc_n = hfar;
  lc_n = std::max(lc_n, CTX::instance()->mesh.lcMin);
  lc_n = std::min(lc_n, CTX::instance()->mesh.lcMax);
  SVector3 t1= SVector3(x-xpk,y-ypk,z-zpk);
  t1.normalize();
  metr = buildMetricTangentToCurve(t1,lc_n,lc_n);
}

void BoundaryLayerField::operator() (AttractorField *cc, double dist,
                                     double x, double y, double z,
                                     SMetric3 &metr, GEntity *ge)
{
  // dist = hwall -> lc = hwall * ratio
  // dist = hwall (1+ratio) -> lc = hwall ratio ^ 2