Skip to content
Snippets Groups Projects
Select Git revision
  • 509e02db33faf9d6cc2a19112382e765ed5b9553
  • master default
  • cgnsUnstructured
  • partitioning
  • poppler
  • HighOrderBLCurving
  • gmsh_3_0_4
  • gmsh_3_0_3
  • gmsh_3_0_2
  • gmsh_3_0_1
  • gmsh_3_0_0
  • gmsh_2_16_0
  • gmsh_2_15_0
  • gmsh_2_14_1
  • gmsh_2_14_0
  • gmsh_2_13_2
  • gmsh_2_13_1
  • gmsh_2_12_0
  • gmsh_2_11_0
  • gmsh_2_10_1
  • gmsh_2_10_0
  • gmsh_2_9_3
  • gmsh_2_9_2
  • gmsh_2_9_1
  • gmsh_2_9_0
  • gmsh_2_8_6
26 results

xmath.cpp

Blame
  • Forked from gmsh / gmsh
    Source project has a limited visibility.
    • Christophe Geuzaine's avatar
      504b5126
      · 504b5126
      Christophe Geuzaine authored
      - removed all the crappy STL code and rewrote it using JF's POLY_rep
        class
      
      - generalized POLY_rep so that we can use the polygonal discretization
        as a surface mesh, and mesh in 3D afterwards. I.e., we can now take
        an input triangulation (a single surface in STL format, multiple
        surfaces in STL format, one or more surfaces defined using the new
        "Discrete Surface" commands), and generate a 3D mesh that uses it. We
        could in theory even mix triangulated and "normal" surfaces in the
        same geometry, but nothing is done at the moment to ensure that the
        mesh at the interfaces would match (if it does, it actually works
        very nicely)
      
      - new STL mesh output format to export a surface mesh as a STL file
      
      - added an option to the GEO output routine to save the surface mesh
        as discrete surfaces associated with the geometrical surfaces
      
      - added STL and Text output formats for post-processing views (the
        text output allows for example to exploit plugin-generated data in
        gnuplot)
      
      - generalized Plugin(Evaluate):
      
        * can loop automatically over all the timestep and/or components
      
        * can do operations using data from an external view
      
           - if the 2 views are based on the same grid, the plugin does the
             evaluation very efficiently
      
           - if the 2 views are based on differenet grids, the plugin
             automatically interpolates the external view data onto the
             grid of the current view
      
      - added new Rand() function in MathEval
      
      - default colormap is now # 2 (the Matlab "Jet" colormap)
      504b5126
      History
      Christophe Geuzaine authored
      - removed all the crappy STL code and rewrote it using JF's POLY_rep
        class
      
      - generalized POLY_rep so that we can use the polygonal discretization
        as a surface mesh, and mesh in 3D afterwards. I.e., we can now take
        an input triangulation (a single surface in STL format, multiple
        surfaces in STL format, one or more surfaces defined using the new
        "Discrete Surface" commands), and generate a 3D mesh that uses it. We
        could in theory even mix triangulated and "normal" surfaces in the
        same geometry, but nothing is done at the moment to ensure that the
        mesh at the interfaces would match (if it does, it actually works
        very nicely)
      
      - new STL mesh output format to export a surface mesh as a STL file
      
      - added an option to the GEO output routine to save the surface mesh
        as discrete surfaces associated with the geometrical surfaces
      
      - added STL and Text output formats for post-processing views (the
        text output allows for example to exploit plugin-generated data in
        gnuplot)
      
      - generalized Plugin(Evaluate):
      
        * can loop automatically over all the timestep and/or components
      
        * can do operations using data from an external view
      
           - if the 2 views are based on the same grid, the plugin does the
             evaluation very efficiently
      
           - if the 2 views are based on differenet grids, the plugin
             automatically interpolates the external view data onto the
             grid of the current view
      
      - added new Rand() function in MathEval
      
      - default colormap is now # 2 (the Matlab "Jet" colormap)
    Field.cpp 62.06 KiB
    // Gmsh - Copyright (C) 1997-2011 C. Geuzaine, J.-F. Remacle
    //
    // See the LICENSE.txt file for license information. Please report all
    // bugs and problems to <gmsh@geuz.org>.
    //
    // Contributor(s):
    //   Jonathan Lambrechts
    //
    
    #include <stdlib.h>
    #include <list>
    #include <math.h>
    #include <fstream>
    #include <string>
    #include <string.h>
    #include <sstream>
    #include "GmshConfig.h"
    #include "Context.h"
    #include "Field.h"
    #include "GeoInterpolation.h"
    #include "GModel.h"
    #include "GmshMessage.h"
    #include "Numeric.h"
    #include "mathEvaluator.h"
    #include "BackgroundMesh.h"
    
    #if defined(HAVE_POST)
    #include "OctreePost.h"
    #include "PViewDataList.h"
    #include "MVertex.h"
    #endif
    
    #if defined(HAVE_ANN)
    #include "ANN/ANN.h"
    #endif
    
    class FieldOptionDouble : public FieldOption
    {
     public:
      double &val;
      FieldOptionType getType(){ return FIELD_OPTION_DOUBLE; }
      FieldOptionDouble(double &_val, std::string _help, bool *_status=0)
        : FieldOption(_help, _status), val(_val){}
      double numericalValue() const { return val; }
      void numericalValue(double v){ modified(); val = v; }
      void getTextRepresentation(std::string &v_str)
      {
        std::ostringstream sstream;
        sstream.precision(16);
        sstream << val;
        v_str = sstream.str();
      }
    };
    
    class FieldOptionInt : public FieldOption
    {
     public:
      int &val;
      FieldOptionType getType(){ return FIELD_OPTION_INT; }
      FieldOptionInt(int &_val, std::string _help, bool *_status=0) 
        : FieldOption(_help, _status), val(_val){}
      double numericalValue() const { return val; }
      void numericalValue(double v){ modified(); val = (int)v; }
      void getTextRepresentation(std::string & v_str)
      {
        std::ostringstream sstream;
        sstream << val;
        v_str = sstream.str();
      }
    };
    
    class FieldOptionList : public FieldOption
    {
     public:
      std::list<int> &val;
      FieldOptionType getType(){ return FIELD_OPTION_LIST; }
      FieldOptionList(std::list<int> &_val, std::string _help, bool *_status=0) 
        : FieldOption(_help, _status), val(_val) {}
      std::list<int> &list(){ modified(); return val; }
      const std::list<int>& list() const { return val; }
      void getTextRepresentation(std::string & v_str)
      {
        std::ostringstream sstream;
        sstream << "{";
        for(std::list<int>::iterator it = val.begin(); it != val.end(); it++) {
          if(it != val.begin())
            sstream << ", ";
          sstream << *it;
        }
        sstream << "}";
        v_str = sstream.str();
      }
    };
    
    class FieldOptionString : public FieldOption
    {
     public:
      std::string & val;
      virtual FieldOptionType getType(){ return FIELD_OPTION_STRING; }
      FieldOptionString(std::string &_val, std::string _help, bool *_status=0)
        : FieldOption(_help, _status), val(_val) {}
      std::string &string() { modified(); return val; }
      const std::string &string() const { return val; }
      void getTextRepresentation(std::string &v_str)
      {
        std::ostringstream sstream;
        sstream << "\"" << val << "\"";
        v_str = sstream.str();
      }
    };
    
    class FieldOptionPath : public FieldOptionString
    {
     public:
      virtual FieldOptionType getType(){ return FIELD_OPTION_PATH; }
      FieldOptionPath(std::string &_val, std::string _help, bool *_status=0)
        : FieldOptionString(_val, _help, _status) {}
    };
    
    class FieldOptionBool : public FieldOption
    {
     public:
      bool & val;
      FieldOptionType getType(){ return FIELD_OPTION_BOOL; }
      FieldOptionBool(bool & _val, std::string _help, bool *_status=0)
        : FieldOption(_help, _status), val(_val) {}
      double numericalValue() const { return val; }
      void numericalValue(double v){ modified(); val = v; }
      void getTextRepresentation(std::string & v_str)
      {
        std::ostringstream sstream;
        sstream << val;
        v_str = sstream.str();
      }
    };
    
    void FieldManager::reset()
    {
      for(std::map<int, Field *>::iterator it = begin(); it != end(); it++) {
        delete it->second;
      }
      clear();
    }
    
    Field *FieldManager::get(int id)
    {
      iterator it = find(id);
      if(it == end()) return NULL;
      return it->second;
    }
    
    Field *FieldManager::newField(int id, std::string type_name)
    {
      if(find(id) != end()) {
        Msg::Error("Field id %i is already defined", id);
        return 0;
      }
      if(map_type_name.find(type_name) == map_type_name.end()) {
        Msg::Error("Unknown field type \"%s\"", type_name.c_str());
        return 0;
      }
      Field *f = (*map_type_name[type_name]) ();
      if(!f)
        return 0;
      f->id = id;
      (*this)[id] = f;
      return f;
    }
    
    int FieldManager::newId()
    {
      int i = 0;
      iterator it = begin();
      while(1) {
        i++;
        while(it != end() && it->first < i)
          it++;
        if(it == end() || it->first != i)
          break;
      }
      return std::max(i, 1);
    }
    
    int FieldManager::maxId()
    {
      if(!empty())
        return rbegin()->first;
      else
        return 0;
    }
    
    void FieldManager::deleteField(int id)
    {
      iterator it = find(id);
      if(it == end()) {
        Msg::Error("Cannot delete field id %i, it does not exist", id);
        return;
      }
      delete it->second;
      erase(it);
    }
    
    // StructuredField
    class StructuredField : public Field
    {
      double o[3], d[3];
      int n[3];
      double *data;
      bool error_status;
      bool text_format;
      std::string file_name;
     public:
      StructuredField()
      {
        options["FileName"] = new FieldOptionPath
          (file_name, "Name of the input file", &update_needed);
        text_format = false;
        options["TextFormat"] = new FieldOptionBool
          (text_format, "True for ASCII input files, false for binary files (4 bite "
           "signed integers for n, double precision floating points for v, D and O)", 
           &update_needed);
        data = 0;
      }
      std::string getDescription()
      {
        return "Linearly interpolate between data provided on a 3D rectangular "
          "structured grid.\n\n"
          "The format of the input file is:\n\n"
          "  Ox Oy Oz \n"
          "  Dx Dy Dz \n"
          "  nx ny nz \n"
          "  v(0,0,0) v(0,0,1) v(0,0,2) ... \n"
          "  v(0,1,0) v(0,1,1) v(0,1,2) ... \n"
          "  v(0,2,0) v(0,2,1) v(0,2,2) ... \n"
          "  ...      ...      ... \n"
          "  v(1,0,0) ...      ... \n\n"
          "where O are the coordinates of the first node, D are the distances "
          "between nodes in each direction, n are the numbers of nodes in each "
          "direction, and v are the values on each node.";
      }
      const char *getName()
      {
        return "Structured";
      }
      virtual ~StructuredField()
      {
        if(data) delete[]data;
      }
      double operator() (double x, double y, double z, GEntity *ge=0)
      {
        if(update_needed) {
          error_status = false;
          try {
            std::ifstream input;
            if(text_format)
              input.open(file_name.c_str());
            else
              input.open(file_name.c_str(),std::ios::binary);
            if(!input.is_open())
              throw(1);
            input.
              exceptions(std::ifstream::eofbit | std::ifstream::failbit | std::
                         ifstream::badbit);
            if(!text_format) {
              input.read((char *)o, 3 * sizeof(double));
              input.read((char *)d, 3 * sizeof(double));
              input.read((char *)n, 3 * sizeof(int));
              int nt = n[0] * n[1] * n[2];
              if(data)
                delete[]data;
              data = new double[nt];
              input.read((char *)data, nt * sizeof(double));
            }
            else {
              input >> o[0] >> o[1] >> o[2] >> d[0] >> d[1] >> d[2] >> n[0] >>
                n[1] >> n[2];
              int nt = n[0] * n[1] * n[2];
              if(data)
                delete[]data;
              data = new double[nt];
              for(int i = 0; i < nt; i++)
                input >> data[i];
            }
            input.close();
          }
          catch(...) {
            error_status = true;
            Msg::Error("Field %i : error reading file %s", this->id, file_name.c_str());
          }
          update_needed = false;
        }
        if(error_status)
          return MAX_LC;
        //tri-linear
        int id[2][3];
        double xi[3];
        double xyz[3] = { x, y, z };
        for(int i = 0; i < 3; i++) {
          id[0][i] = (int)floor((xyz[i] - o[i]) / d[i]);
          id[1][i] = id[0][i] + 1;
          id[0][i] = std::max(std::min(id[0][i], n[i] - 1), 0);
          id[1][i] = std::max(std::min(id[1][i], n[i] - 1), 0);
          xi[i] = (xyz[i] - (o[i] + id[0][i] * d[i])) / d[i];
          xi[i] = std::max(std::min(xi[i], 1.), 0.);
        }
        double v = 0;
        for(int i = 0; i < 2; i++)
          for(int j = 0; j < 2; j++)
            for(int k = 0; k < 2; k++) {
              v += data[id[i][0] * n[1] * n[2] + id[j][1] * n[2] + id[k][2]]
                * (i * xi[0] + (1 - i) * (1 - xi[0]))
                * (j * xi[1] + (1 - j) * (1 - xi[1]))
                * (k * xi[2] + (1 - k) * (1 - xi[2]));
            }
        return v;
      }
    };
    
    class UTMField : public Field
    {
      int iField, zone;
      double a, b, n, n2, n3, n4, n5, e, e2, e1, e12, e13, e14, J1, J2, J3, J4,
        Ap, Bp, Cp, Dp, Ep, e4, e6, ep, ep2, ep4, k0, mu_fact;
     public:
      std::string getDescription()
      {
        return "Evaluate Field[IField] in Universal Transverse Mercator coordinates.\n\n"
          "The formulas for the coordinates transformation are taken from:\n\n"
          "  http://www.uwgb.edu/dutchs/UsefulData/UTMFormulas.HTM";
      }
      UTMField()
      {
        iField = 1;
        zone = 0;
        options["IField"] = new FieldOptionInt
          (iField, "Index of the field to evaluate");
        options["Zone"] = new FieldOptionInt
          (zone, "Zone of the UTM projection");
        a = 6378137; // Equatorial Radius
        b = 6356752.3142; // Rayon Polar Radius
        // see http://www.uwgb.edu/dutchs/UsefulData/UTMFormulas.HTM
        n = (a - b) / (a + b);
        n2 = n * n;
        n3 = n * n * n;
        n4 = n * n * n * n;
        n5 = n * n * n * n * n;
        e = sqrt(1 - b * b / a / a);
        e2 = e * e;
        e1 = (1 - sqrt(1 - e2)) / (1 + sqrt(1 - e2));
        e12 = e1 * e1;
        e13 = e1 * e1 * e1;
        e14 = e1 * e1 * e1 * e1;
        J1 = (3 * e1 / 2 - 27 * e13 / 32);
        J2 = (21 * e12 / 16 - 55 * e14 / 32);
        J3 = 151 * e13 / 96;
        J4 = 1097 * e14 / 512;
        Ap = a * (1 - n + (5. / 4.) * (n2 - n3) + (81. / 64.) * (n4 - n5));
        Bp = -3 * a * n / 2 * (1 - n + (7. / 8.) * (n2 - n3) +
                               (55. / 64.) * (n4 - n5));
        Cp = 14 * a * n2 / 16 * (1 - n + (3. / 4) * (n2 - n3));
        Dp = -35 * a * n3 / 48 * (1 - n + 11. / 16. * (n2 - n3));
        Ep = +315 * a * n4 / 51 * (1 - n);
        e4 = e2 * e2;
        e6 = e2 * e2 * e2;
        ep = e * a / b;
        ep2 = ep * ep;
        ep4 = ep2 * ep2;
        k0 = 0.9996;
        mu_fact = 1 / (k0 * a * (1 - e2 / 4 - 3 * e4 / 64 - 5 * e6 / 256));
      }
      const char *getName()
      {
        return "UTM";
      }
      double operator() (double x, double y, double z, GEntity *ge=0)
      {
        double r = sqrt(x * x + y * y + z * z);
        double lon = atan2(y, x);
        double lat = asin(z / r);
        double meridionalarc = Ap * lat + Bp * sin(2 * lat)
          + Cp * sin(4 * lat) + Dp * sin(6 * lat) + Ep;
        double slat = sin(lat);
        double clat = cos(lat);
        double slat2 = slat * slat;
        double clat2 = clat * clat;
        double clat3 = clat2 * clat;
        double clat4 = clat3 * clat;
        double tlat2 = slat2 / clat2;
        double nu = a / sqrt(1 - e * e * slat2);
        double p = lon - ((zone - 0.5) / 30 - 1) * M_PI;
        double p2 = p * p;
        double p3 = p * p2;
        double p4 = p2 * p2;
        double utm_x =
          k0 * nu * clat * p + (k0 * nu * clat3 / 6) * (1 - tlat2 +
                                                        ep2 * clat2) * p3 + 5e5;
        double utm_y =
          meridionalarc * k0 + k0 * nu * slat * clat / 2 * p2 +
          k0 * nu * slat * clat3 / 24 * (5 - tlat2 + 9 * ep2 * clat2 +
                                         4 * ep4 * clat4) * p4;
        Field *field = GModel::current()->getFields()->get(iField);
        if(!field || iField == id) return MAX_LC;
        return (*field)(utm_x, utm_y, 0);
      }
    };
    
    class LonLatField : public Field
    {
      int iField, fromStereo;
      double stereoRadius;
     public:
      std::string getDescription()
      {
        return "Evaluate Field[IField] in geographic coordinates (longitude, latitude):\n\n"
          "  F = Field[IField](atan(y/x), asin(z/sqrt(x^2+y^2+z^2))";
      }
      LonLatField()
      {
        iField = 1;
        options["IField"] = new FieldOptionInt
          (iField, "Index of the field to evaluate.");
        fromStereo = 0;
        stereoRadius = 6371e3;
        
        options["FromStereo"] = new FieldOptionInt
          (fromStereo, "if = 1, the mesh is in stereographic coordinates. "
           "xi = 2Rx/(R+z),  eta = 2Ry/(R+z)");
        options["RadiusStereo"] = new FieldOptionDouble
          (stereoRadius, "radius of the sphere of the stereograpic coordinates");
      }
      const char *getName()
      {
        return "LonLat";
      }
      double operator() (double x, double y, double z, GEntity *ge=0)
      {
        Field *field = GModel::current()->getFields()->get(iField);
        if(!field || iField == id) return MAX_LC;
        if (fromStereo == 1) {
          double xi = x;
          double eta = y;
          double r2 = stereoRadius * stereoRadius;
          x = 4 * r2 * xi / ( 4 * r2 + xi * xi + eta * eta);
          y = 4 * r2 * eta / ( 4 * r2 + xi * xi + eta * eta);
          z = stereoRadius * (4 * r2 - eta * eta - xi * xi) / ( 4 * r2 + xi * xi + eta * eta);
        }
        return (*field)(atan2(y, x), asin(z / stereoRadius), 0);
      }
    };
    
    class BoxField : public Field
    {
      double v_in, v_out, x_min, x_max, y_min, y_max, z_min, z_max;
     public:
      std::string getDescription()
      {
        return "The value of this field is VIn inside the box, VOut outside the box. "
          "The box is given by\n\n"
          "  Xmin <= x <= XMax &&\n"
          "  YMin <= y <= YMax &&\n"
          "  ZMin <= z <= ZMax";
      }
      BoxField()
      {
        v_in = v_out = x_min = x_max = y_min = y_max = z_min = z_max = 0;
        options["VIn"] = new FieldOptionDouble
          (v_in, "Value inside the box");
        options["VOut"] = new FieldOptionDouble
          (v_out, "Value outside the box");
        options["XMin"] = new FieldOptionDouble
          (x_min, "Minimum X coordinate of the box");
        options["XMax"] = new FieldOptionDouble
          (x_max, "Maximum X coordinate of the box");
        options["YMin"] = new FieldOptionDouble
          (y_min, "Minimum Y coordinate of the box");
        options["YMax"] = new FieldOptionDouble
          (y_max, "Maximum Y coordinate of the box");
        options["ZMin"] = new FieldOptionDouble
          (z_min, "Minimum Z coordinate of the box");
        options["ZMax"] = new FieldOptionDouble
          (z_max, "Maximum Z coordinate of the box");    
      }
      const char *getName()
      {
        return "Box";
      }
      double operator() (double x, double y, double z, GEntity *ge=0)
      {    
        return (x <= x_max && x >= x_min && y <= y_max && y >= y_min && z <= z_max
                && z >= z_min) ? v_in : v_out;
      }
    };
    
    class CylinderField : public Field
    {
      double v_in, v_out;
      double xc,yc,zc;
      double xa,ya,za;
      double R;
      
     public:
      std::string getDescription()
      {
        return "The value of this field is VIn inside a frustrated cylinder, VOut outside. "
          "The cylinder is given by\n\n"
          "  ||dX||^2 < R^2 &&\n"
          "  (X-X0).A < ||A||^2\n"
          "  dX = (X - X0) - ((X - X0).A)/(||A||^2) . A";
      }
      CylinderField()
      {
        v_in = v_out = xc = yc = zc = xa = ya = R = 0;
        za = 1.;
        
        options["VIn"] = new FieldOptionDouble
          (v_in, "Value inside the cylinder");
        options["VOut"] = new FieldOptionDouble
          (v_out, "Value outside the cylinder");
    
        options["XCenter"] = new FieldOptionDouble
          (xc, "X coordinate of the cylinder center");
        options["YCenter"] = new FieldOptionDouble
          (yc, "Y coordinate of the cylinder center");
        options["ZCenter"] = new FieldOptionDouble
          (zc, "Z coordinate of the cylinder center");
    
        
        options["XAxis"] = new FieldOptionDouble
          (xa, "X component of the cylinder axis");
        options["YAxis"] = new FieldOptionDouble
          (ya, "Y component of the cylinder axis");
        options["ZAxis"] = new FieldOptionDouble
          (za, "Z component of the cylinder axis");
    
        options["Radius"] = new FieldOptionDouble
          (R,"Radius");    
        
      }
      const char *getName()
      {
        return "Cylinder";
      }
      double operator() (double x, double y, double z, GEntity *ge=0)
      {    
        double dx = x-xc;
        double dy = y-yc;
        double dz = z-zc;
        
        double adx = (xa * dx + ya * dy + za * dz)/(xa*xa + ya*ya + za*za);
        
        dx -= adx * xa;
        dy -= adx * ya;
        dz -= adx * za;
        
        return ((dx*dx + dy*dy + dz*dz < R*R) && fabs(adx) < 1) ? v_in : v_out;
      }
    };
    
    class ThresholdField : public Field
    {
     protected :
      int iField;
      double dmin, dmax, lcmin, lcmax;
      bool sigmoid, stopAtDistMax;
     public:
      virtual const char *getName()
      {
        return "Threshold";
      }
      virtual std::string getDescription()
      {
        return "F = LCMin if Field[IField] <= DistMin,\n"
          "F = LCMax if Field[IField] >= DistMax,\n"
          "F = interpolation between LcMin and LcMax if DistMin < Field[IField] < DistMax";
      }
      ThresholdField()
      {
        iField = 0;
        dmin = 1;
        dmax = 10;
        lcmin = 0.1;
        lcmax = 1;
        sigmoid = false;
        stopAtDistMax = false;
        options["IField"] = new FieldOptionInt
          (iField, "Index of the field to evaluate");
        options["DistMin"] = new FieldOptionDouble
          (dmin, "Distance from entity up to which element size will be LcMin");
        options["DistMax"] = new FieldOptionDouble
          (dmax, "Distance from entity after which element size will be LcMax");
        options["LcMin"] = new FieldOptionDouble
          (lcmin, "Element size inside DistMin");
        options["LcMax"] = new FieldOptionDouble
          (lcmax, "Element size outside DistMax");
        options["Sigmoid"] = new FieldOptionBool
          (sigmoid, "True to interpolate between LcMin and LcMax using a sigmoid, "
           "false to interpolate linearly");
        options["StopAtDistMax"] = new FieldOptionBool
          (stopAtDistMax, "True to not impose element size outside DistMax (i.e., "
           "F = a very big value if Field[IField] > DistMax)");
      }
      double operator() (double x, double y, double z, GEntity *ge=0)
      {
        Field *field = GModel::current()->getFields()->get(iField);
        if(!field || iField == id) return MAX_LC;
        double r = ((*field) (x, y, z) - dmin) / (dmax - dmin);
        r = std::max(std::min(r, 1.), 0.);
        double lc;
        if(stopAtDistMax && r >= 1.){
          lc = MAX_LC;
        }
        else if(sigmoid){
          double s = exp(12. * r - 6.) / (1. + exp(12. * r - 6.));
          lc = lcmin * (1. - s) + lcmax * s;
        }
        else{ // linear
          lc = lcmin * (1 - r) + lcmax * r;
        }
        return lc;
      }
    };
    
    class GradientField : public Field
    {
      int iField, kind;
      double delta;
     public:
      const char *getName()
      {
        return "Gradient";
      }
      std::string getDescription()
      {
        return "Compute the finite difference gradient of Field[IField]:\n\n"
          "  F = (Field[IField](X + Delta/2) - "
          "       Field[IField](X - Delta/2)) / Delta";
      }
      GradientField() : iField(0), kind(3), delta(CTX::instance()->lc / 1e4)
      {
        iField = 1;
        kind = 0;
        delta = 0.;
        options["IField"] = new FieldOptionInt
          (iField, "Field index");
        options["Kind"] = new FieldOptionInt
          (kind, "Component of the gradient to evaluate: 0 for X, 1 for Y, 2 for Z, "
           "3 for the norm");
        options["Delta"] = new FieldOptionDouble
          (delta, "Finite difference step");
      }
      double operator() (double x, double y, double z, GEntity *ge=0)
      {
        Field *field = GModel::current()->getFields()->get(iField);
        if(!field || iField == id) return MAX_LC;
        double gx, gy, gz;
        switch (kind) {
        case 0:    /* x */
          return ((*field) (x + delta / 2, y, z) -
                  (*field) (x - delta / 2, y, z)) / delta;
        case 1:    /* y */
          return ((*field) (x, y + delta / 2, z) -
                  (*field) (x, y - delta / 2, z)) / delta;
        case 2:    /* z */
          return ((*field) (x, y, z + delta / 2) -
                  (*field) (x, y, z - delta / 2)) / delta;
        case 3:    /* norm */
          gx =
            ((*field) (x + delta / 2, y, z) -
             (*field) (x - delta / 2, y, z)) / delta;
          gy =
            ((*field) (x, y + delta / 2, z) -
             (*field) (x, y - delta / 2, z)) / delta;
          gz =
            ((*field) (x, y, z + delta / 2) -
             (*field) (x, y, z - delta / 2)) / delta;
          return sqrt(gx * gx + gy * gy + gz * gz);
        default:
          Msg::Error("Field %i : Unknown kind (%i) of gradient", this->id,
                     kind);
          return MAX_LC;
        }
      }
    };
    
    class CurvatureField : public Field
    {
      int iField;
      double delta;
     public:
      const char *getName()
      {
        return "Curvature";
      }
      std::string getDescription()
      {
        return "Compute the curvature of Field[IField]:\n\n"
          "  F = div(norm(grad(Field[IField])))";
      }
      CurvatureField() : iField(0), delta(CTX::instance()->lc / 1e4)
      {
        iField = 1;
        delta = 0.;
        options["IField"] = new FieldOptionInt
          (iField, "Field index");
        options["Delta"] = new FieldOptionDouble
          (delta, "Step of the finite differences");
      }
      void grad_norm(Field &f, double x, double y, double z, double *g)
      {
        g[0] = f(x + delta / 2, y, z) - f(x - delta / 2, y, z);
        g[1] = f(x, y + delta / 2, z) - f(x, y - delta / 2, z);
        g[2] = f(x, y, z + delta / 2) - f(x, y, z - delta / 2);
        double n=sqrt(g[0] * g[0] + g[1] * g[1] + g[2] * g[2]);
        g[0] /= n;
        g[1] /= n;
        g[2] /= n;
      }
      double operator() (double x, double y, double z, GEntity *ge=0)
      {
        Field *field = GModel::current()->getFields()->get(iField);
        if(!field || iField == id) return MAX_LC;
        double grad[6][3];
        grad_norm(*field, x + delta / 2, y, z, grad[0]);
        grad_norm(*field, x - delta / 2, y, z, grad[1]);
        grad_norm(*field, x, y + delta / 2, z, grad[2]);
        grad_norm(*field, x, y - delta / 2, z, grad[3]);
        grad_norm(*field, x, y, z + delta / 2, grad[4]);
        grad_norm(*field, x, y, z - delta / 2, grad[5]);
        return (grad[0][0] - grad[1][0] + grad[2][1] - 
                grad[3][1] + grad[4][2] - grad[5][2]) / delta;
      }
    };
    
    class MaxEigenHessianField : public Field
    {
      int iField;
      double delta;
     public:
      const char *getName()
      {
        return "MaxEigenHessian";
      }
      std::string getDescription()
      {
        return "Compute the maximum eigenvalue of the Hessian matrix of "
          "Field[IField], with the gradients evaluated by finite differences:\n\n"
          "  F = max(eig(grad(grad(Field[IField]))))";
      }
      MaxEigenHessianField() : iField(0), delta(CTX::instance()->lc / 1e4)
      {
        iField = 1;
        delta = 0.;
        options["IField"] = new FieldOptionInt
          (iField, "Field index");
        options["Delta"] = new FieldOptionDouble
          (delta, "Step used for the finite differences");
      }
      double operator() (double x, double y, double z, GEntity *ge=0)
      {
        Field *field = GModel::current()->getFields()->get(iField);
        if(!field || iField == id) return MAX_LC;
        double mat[3][3], eig[3];
        mat[1][0] = mat[0][1] = (*field) (x + delta / 2 , y + delta / 2, z)
          + (*field) (x - delta / 2 , y - delta / 2, z)
          - (*field) (x - delta / 2 , y + delta / 2, z)
          - (*field) (x + delta / 2 , y - delta / 2, z);
        mat[2][0] = mat[0][2] = (*field) (x + delta/2 , y, z + delta / 2)
          + (*field) (x - delta / 2 , y, z - delta / 2)
          - (*field) (x - delta / 2 , y, z + delta / 2)
          - (*field) (x + delta / 2 , y, z - delta / 2);
        mat[2][1] = mat[1][2] = (*field) (x, y + delta/2 , z + delta / 2)
          + (*field) (x, y - delta / 2 , z - delta / 2)
          - (*field) (x, y - delta / 2 , z + delta / 2)
          - (*field) (x, y + delta / 2 , z - delta / 2);
        double f = (*field)(x, y, z);
        mat[0][0] = (*field)(x + delta, y, z) + (*field)(x - delta, y, z) - 2 * f;
        mat[1][1] = (*field)(x, y + delta, z) + (*field)(x, y - delta, z) - 2 * f;
        mat[2][2] = (*field)(x, y, z + delta) + (*field)(x, y, z - delta) - 2 * f;
        eigenvalue(mat, eig);
        return eig[0] / (delta * delta);
      }
    };
    
    class LaplacianField : public Field
    {
      int iField;
      double delta;
     public:
      const char *getName()
      {
        return "Laplacian";
      }
      std::string getDescription()
      {
        return "Compute finite difference the Laplacian of Field[IField]:\n\n"
          "  F = G(x+d,y,z) + G(x-d,y,z) +\n"
          "      G(x,y+d,z) + G(x,y-d,z) +\n"
          "      G(x,y,z+d) + G(x,y,z-d) - 6 * G(x,y,z),\n\n"
          "where G=Field[IField] and d=Delta";
      }
      LaplacianField() : iField(0), delta(CTX::instance()->lc / 1e4)
      {
        iField = 1;
        delta = 0.1;
        options["IField"] = new FieldOptionInt
          (iField, "Field index");
        options["Delta"] = new FieldOptionDouble
          (delta, "Finite difference step");
      }
      double operator() (double x, double y, double z, GEntity *ge=0)
      {
        Field *field = GModel::current()->getFields()->get(iField);
        if(!field || iField == id) return MAX_LC;
        return ((*field) (x + delta , y, z)+ (*field) (x - delta , y, z)
                +(*field) (x, y + delta , z)+ (*field) (x, y - delta , z)
                +(*field) (x, y, z + delta )+ (*field) (x, y, z - delta )
                -6* (*field) (x , y, z)) / (delta*delta);
      }
    };
    
    class MeanField : public Field
    {
      int iField;
      double delta;
     public:
      const char *getName()
      {
        return "Mean";
      }
      std::string getDescription()
      {
        return "Simple smoother:\n\n"
          "  F = (G(x+delta,y,z) + G(x-delta,y,z) +\n"
          "       G(x,y+delta,z) + G(x,y-delta,z) +\n"
          "       G(x,y,z+delta) + G(x,y,z-delta) +\n"
          "       G(x,y,z)) / 7,\n\n"
          "where G=Field[IField]";
      }
      MeanField() : iField(0), delta(CTX::instance()->lc / 1e4)
      {
        options["IField"] = new FieldOptionInt
          (iField, "Field index");
        options["Delta"] = new FieldOptionDouble
          (delta, "Distance used to compute the mean value");
      }
      double operator() (double x, double y, double z, GEntity *ge=0)
      {
        Field *field = GModel::current()->getFields()->get(iField);
        if(!field || iField == id) return MAX_LC;
        return ((*field) (x + delta , y, z) + (*field) (x - delta, y, z)
                + (*field) (x, y + delta, z) + (*field) (x, y - delta, z)
                + (*field) (x, y, z + delta) + (*field) (x, y, z - delta)
                + (*field) (x, y, z)) / 7;
      }
    };
    
    class MathEvalExpression
    {
     private:
      mathEvaluator *_f;
      std::set<int> _fields;
     public:
      MathEvalExpression() : _f(0) {}
      ~MathEvalExpression(){ if(_f) delete _f; }
      bool set_function(const std::string &f)
      {
        // get id numbers of fields appearing in the function
        _fields.clear();
        unsigned int i = 0;
        while(i < f.size()){
          unsigned int j = 0;
          if(f[i] == 'F'){
            std::string id("");
            while(i + 1 + j < f.size() && f[i + 1 + j] >= '0' && f[i + 1 + j] <= '9'){
              id += f[i + 1 + j];
              j++;
            }
            _fields.insert(atoi(id.c_str()));
          }
          i += j + 1;
        }
        std::vector<std::string> expressions(1), variables(3 + _fields.size());
        expressions[0] = f;
        variables[0] = "x"; 
        variables[1] = "y"; 
        variables[2] = "z";
        i = 3;
        for(std::set<int>::iterator it = _fields.begin(); it != _fields.end(); it++){
          std::ostringstream sstream;
          sstream << "F" << *it;
          variables[i++] = sstream.str();
        }
        if(_f) delete _f;
        _f = new mathEvaluator(expressions, variables);
        if(expressions.empty()) {
          delete _f;
          _f = 0;
          return false;
        }
        return true;
      }
      double evaluate(double x, double y, double z)
      {
        if(!_f) return MAX_LC;
        std::vector<double> values(3 + _fields.size()), res(1);
        values[0] = x;
        values[1] = y;
        values[2] = z;
        int i = 3;
        for(std::set<int>::iterator it = _fields.begin(); it != _fields.end(); it++){
          Field *field = GModel::current()->getFields()->get(*it);
          values[i++] = field ? (*field)(x, y, z) : MAX_LC;
        }
        if(_f->eval(values, res)) 
          return res[0];
        else
          return MAX_LC;
      }
    };
    
    class MathEvalExpressionAniso
    {
     private:
      mathEvaluator *_f[6];
      std::set<int> _fields[6];
     public:
      MathEvalExpressionAniso() {
        for (int i=0;i<6;i++)_f[i]=0; 
      }
      ~MathEvalExpressionAniso(){ 
        for (int i=0;i<6;i++)if(_f[i]) delete _f[i]; 
      }
      bool set_function(int iFunction, const std::string &f)
      {
        // get id numbers of fields appearing in the function
        _fields[iFunction].clear();
        unsigned int i = 0;
        while(i < f.size()){
          unsigned int j = 0;
          if(f[i] == 'F'){
            std::string id("");
            while(i + 1 + j < f.size() && f[i + 1 + j] >= '0' && f[i + 1 + j] <= '9'){
              id += f[i + 1 + j];
              j++;
            }
            _fields[iFunction].insert(atoi(id.c_str()));
          }
          i += j + 1;
        }
        std::vector<std::string> expressions(1), variables(3 + _fields[iFunction].size());
        expressions[0] = f;
        variables[0] = "x"; 
        variables[1] = "y"; 
        variables[2] = "z";
        i = 3;
        for(std::set<int>::iterator it = _fields[iFunction].begin(); 
            it != _fields[iFunction].end(); it++){
          std::ostringstream sstream;
          sstream << "F" << *it;
          variables[i++] = sstream.str();
        }
        if(_f[iFunction]) delete _f[iFunction];
        _f[iFunction] = new mathEvaluator(expressions, variables);
        if(expressions.empty()) {
          delete _f[iFunction];
          _f[iFunction] = 0;
          return false;
        }
        return true;
      }
      void evaluate (double x, double y, double z, SMetric3 &metr)
      {
        const int index[6][2] = {{0,0},{1,1},{2,2},{0,1},{0,2},{1,2}};
        for (int iFunction = 0;iFunction<6;iFunction++){
          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:
      MathEvalField()
      {
        options["F"] = new FieldOptionString
          (f, "Mathematical function to evaluate.", &update_needed);
        f = "F2 + Sin(z)";
      }
      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;
     public:
      int view_index;
      bool crop_negative_values;
      double operator() (double x, double y, double z, GEntity *ge=0)
      {
        // we should maybe test the unique view num instead, but that
        // would be slower
        if(view_index < 0 || view_index >= (int)PView::list.size())
          return MAX_LC;
        PView *v = PView::list[view_index];
        if(v->getData()->hasModel(GModel::current())){
          Msg::Error("Cannot use view based on current model for background mesh");
          Msg::Error("Use a list-based view (.pos file) instead?");
          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 element found containing point (%g,%g,%g)", x, y, z);
        if(l <= 0 && crop_negative_values) return MAX_LC;
        return l;
      }
      const char *getName()
      {
        return "PostView";
      }
      std::string getDescription()
      {
        return "Evaluate the post processing view IView.";
      }
      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;
      }
    };
    #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(v,ff);
          }
        }
        metr = v;
      }
      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 "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;
          }
          int totpoints = nodes_id.size() + (n_nodes_by_edge-2) * edges_id.size() + 
            n_nodes_by_edge * n_nodes_by_edge * faces_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) {
            Vertex *v = FindPoint(*it);
            if(v) {
              getCoord(v->Pos.X, v->Pos.Y, v->Pos.Z, zeronodes[k][0],
                       zeronodes[k][1], zeronodes[k][2]);
    	  _infos[k++] = AttractorInfo(*it,0,0,0);
            }
            else {
              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) {
            Curve *c = FindCurve(*it);
    	GEdge *e = GModel::current()->getEdgeByTag(*it);	
            if(c && !e) {
              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);
                getCoord(V.Pos.X, V.Pos.Y, V.Pos.Z, zeronodes[k][0],
                         zeronodes[k][1], zeronodes[k][2]);
    	    _infos[k++] = AttractorInfo(*it,1,u,0);
              }
            }
            else {
              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);
                  getCoord(gp.x(), gp.y(), gp.z(), zeronodes[k][0],
                           zeronodes[k][1], zeronodes[k][2], e);
    	      _infos[k++] = AttractorInfo(*it,1,u,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.)
          for(std::list<int>::iterator it = faces_id.begin();
              it != faces_id.end(); ++it) {
            Surface *s = FindSurface(*it);
            if(s) {
              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);
                  Vertex V = InterpolateSurface(s, u, v, 0, 0);
                  getCoord(V.Pos.X, V.Pos.Y, V.Pos.Z, zeronodes[k][0],
                           zeronodes[k][1], zeronodes[k][2]);
    	      _infos[k++] = AttractorInfo(*it,2,u,v);
                }
              }
            }
            else {
              GFace *f = GModel::current()->getFaceByTag(*it);
              if(f) {
                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 = ge->parBounds(0);
                    Range<double> b2 = ge->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);
                  }
                }
              }
            }
          }
          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;    
      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["IField"] = new FieldOptionInt
      //    (iField, "Index of the field that contains the distance function");
      options["hwall_n"] = new FieldOptionDouble
        (hwall_n, "Mesh Size Normal to the The Wall");
      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");
    }
    
    double BoundaryLayerField :: operator() (double x, double y, double z, GEntity *ge)
    {
      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);
      //  printf("coucou\n");
      return std::min (hfar,lc);
    }
    
    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
      // dist = hwall (1+ratio+ratio^2) -> lc = hwall ratio ^ 3
      // dist = hwall (1+ratio+ratio^2+...+ratio^{m-1}) = (ratio^{m} - 1)/(ratio-1) 
      // -> lc = hwall ratio ^ m
      // -> find m
      // dist/hwall = (ratio^{m} - 1)/(ratio-1)
      // (dist/hwall)*(ratio-1) + 1 = ratio^{m}
      // lc =  dist*(ratio-1) + hwall 
    
      const double ll1   = dist*(ratio-1) + hwall_n;
      double lc_n  = std::min(ll1,hfar);
      const double ll2   = dist*(ratio-1) + hwall_t;
      double lc_t  = std::min(lc_n*CTX::instance()->mesh.anisoMax, std::min(ll2,hfar));
      
      SVector3 t1,t2,t3;
      double L1,L2,L3;
      
      std::pair<AttractorInfo,SPoint3> pp = cc->getAttractorInfo();
      if (pp.first.dim ==0){
        L1 = lc_n;
        L2 = lc_n;
        L3 = lc_n;
        GVertex *v = GModel::current()->getVertexByTag(pp.first.ent);
        t1 = SVector3(v->x() -x,v->y() -y,v->z() -z);
        t1.normalize();
        //      printf("%p %g %g %g\n",v,t1.x(),t1.y(),t1.z());
        if (t1.norm() == 0)t1 = SVector3(1,0,0);
        if (fabs(t1.x()) < fabs(t1.y()) && fabs(t1.x()) < fabs(t1.z()))
          t2 = SVector3(1,0,0);
        else if (fabs(t1.y()) < fabs(t1.x()) && fabs(t1.y()) < fabs(t1.z()))
          t2 = SVector3(0,1,0);
        else
          t2 = SVector3(0,0,1);
        t3 = crossprod(t1,t2);
        t3.normalize();
        t2 = crossprod(t3,t1);
        //      printf("t1 %p %g %g %g\n",v,t1.x(),t1.y(),t1.z());
        //      printf("t2 %p %g %g %g\n",v,t2.x(),t2.y(),t2.z());
        //      printf("t3 %p %g %g %g\n",v,t3.x(),t3.y(),t3.z());
      }
      else if (pp.first.dim ==1){
        GEdge *e = GModel::current()->getEdgeByTag(pp.first.ent);
        
        // the tangent size at this point is the size of the
        // 1D mesh at this point !
        // hack : use curvature
        /*
          if(CTX::instance()->mesh.lcFromCurvature){
          double Crv = e->curvature(pp.first.u);
          double lc = Crv > 0 ? 2 * M_PI / Crv / CTX::instance()->mesh.minCircPoints : 1.e-22;
          const double ll2b   = dist*(ratio-1) + lc;
          double lc_tb  = std::min(lc_n*CTX::instance()->mesh.anisoMax, std::min(ll2b,hfar));
          lc_t = std::min(lc_t,lc_tb);
          }
          */      
        
        if (dist < hwall_t){
          L1 = lc_t;
          L2 = lc_n;
          L3 = lc_n;
          t1 = e->firstDer(pp.first.u);
          t1.normalize();
          if (fabs(t1.x()) < fabs(t1.y()) && fabs(t1.x()) < fabs(t1.z()))
    	t2 = SVector3(1,0,0);
          else if (fabs(t1.y()) < fabs(t1.x()) && fabs(t1.y()) < fabs(t1.z()))
    	t2 = SVector3(0,1,0);
          else
    	t2 = SVector3(0,0,1);
          t3 = crossprod(t1,t2);
          t3.normalize();
          t2 = crossprod(t3,t1);
          //	  printf("hfar = %g lc = %g dir %g %g \n",hfar,lc,t1.x(),t1.y());
        }
        else {
          L1 = lc_t;
          L2 = lc_n;
          L3 = lc_t;
          GPoint p = e->point(pp.first.u);
          t2 = SVector3(p.x() -x,p.y() -y,p.z() -z);
          if (fabs(t2.x()) < fabs(t2.y()) && fabs(t2.x()) < fabs(t2.z()))
    	  t1 = SVector3(1,0,0);
          else if (fabs(t2.y()) < fabs(t2.x()) && fabs(t2.y()) < fabs(t2.z()))
    	t1 = SVector3(0,1,0);
          else
    	t1 = SVector3(0,0,1);
          t2.normalize();
          t3 = crossprod(t1,t2);
          t3.normalize();
          t1 = crossprod(t3,t2);	  
        }
      }
      else {
        GFace *gf = GModel::current()->getFaceByTag(pp.first.ent);
        
        if (dist < hwall_t){
          L1 = lc_n;
          L2 = lc_t;
          L3 = lc_t;
          t1 = gf->normal(SPoint2(pp.first.u,pp.first.v));
          t1.normalize();
          if (fabs(t1.x()) < fabs(t1.y()) && fabs(t1.x()) < fabs(t1.z()))
    	t2 = SVector3(1,0,0);
          else if (fabs(t1.y()) < fabs(t1.x()) && fabs(t1.y()) < fabs(t1.z()))
    	t2 = SVector3(0,1,0);
          else
    	t2 = SVector3(0,0,1);
          t3 = crossprod(t1,t2);
          t3.normalize();
          t2 = crossprod(t3,t1);
          //	  printf("hfar = %g lc = %g dir %g %g \n",hfar,lc,t1.x(),t1.y());
        }
        else {
          L1 = lc_t;
          L2 = lc_n;
          L3 = lc_t;
          GPoint p = gf->point(SPoint2(pp.first.u,pp.first.v));
          t2 = SVector3(p.x() -x,p.y() -y,p.z() -z);
          if (fabs(t2.x()) < fabs(t2.y()) && fabs(t2.x()) < fabs(t2.z()))
    	t1 = SVector3(1,0,0);
          else if (fabs(t2.y()) < fabs(t2.x()) && fabs(t2.y()) < fabs(t2.z()))
    	t1 = SVector3(0,1,0);
          else
    	t1 = SVector3(0,0,1);
          t2.normalize();
          t3 = crossprod(t1,t2);
          t3.normalize();
          t1 = crossprod(t3,t2);	  
        }	
      }
      metr  = SMetric3(1./(L1*L1), 1./(L2*L2), 1./(L3*L3), t1, t2, t3);
    }
    
    void BoundaryLayerField :: operator() (double x, double y, double z, 
                                           SMetric3 &metr, 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,300));
        }
        update_needed = false;
      }
    
      current_distance = 1.e22;
      current_closest = 0;
      SMetric3 v (1./MAX_LC);
      std::vector<SMetric3> hop;
      for (std::list<AttractorField*>::iterator it = _att_fields.begin();
           it != _att_fields.end(); ++it){
        double cdist = (*(*it)) (x, y, z);
        SPoint3 CLOSEST= (*it)->getAttractorInfo().second;
    
        SMetric3 localMetric;
        (*this)(*it, cdist,x, y, z, localMetric, ge);      
        hop.push_back(localMetric);
        //v = intersection(v,localMetric);
        if (cdist < current_distance){
          current_distance = cdist;
          current_closest = *it;
          v = localMetric;
          _closest_point = CLOSEST;
        }
      }
      //  for (int i=0;i<hop.size();i++)v = intersection_conserveM1(v,hop[i]);
      metr = v;
    }
    #endif
    
    template<class F> class FieldFactoryT : public FieldFactory {
     public:
      Field * operator()() { return new F; }
    };
    
    template<class F> Field *field_factory()
    {
      return new F();
    };
    
    FieldManager::FieldManager()
    {
      map_type_name["Structured"] = new FieldFactoryT<StructuredField>();
      map_type_name["Threshold"] = new FieldFactoryT<ThresholdField>();
    #if defined(HAVE_ANN)
      map_type_name["BoundaryLayer"] = new FieldFactoryT<BoundaryLayerField>();
    #endif
      map_type_name["Box"] = new FieldFactoryT<BoxField>();
      map_type_name["Cylinder"] = new FieldFactoryT<CylinderField>();
      map_type_name["LonLat"] = new FieldFactoryT<LonLatField>();
    #if defined(HAVE_POST)
      map_type_name["PostView"] = new FieldFactoryT<PostViewField>();
    #endif
      map_type_name["Gradient"] = new FieldFactoryT<GradientField>();
      map_type_name["Restrict"] = new FieldFactoryT<RestrictField>();
      map_type_name["Min"] = new FieldFactoryT<MinField>();
      map_type_name["MinAniso"] = new FieldFactoryT<MinAnisoField>();
      map_type_name["Max"] = new FieldFactoryT<MaxField>();
      map_type_name["UTM"] = new FieldFactoryT<UTMField>();
      map_type_name["Laplacian"] = new FieldFactoryT<LaplacianField>();
      map_type_name["Mean"] = new FieldFactoryT<MeanField>();
      map_type_name["Curvature"] = new FieldFactoryT<CurvatureField>();
      map_type_name["Param"] = new FieldFactoryT<ParametricField>();
      map_type_name["MathEval"] = new FieldFactoryT<MathEvalField>();
      map_type_name["MathEvalAniso"] = new FieldFactoryT<MathEvalFieldAniso>();
    #if defined(HAVE_ANN)
      map_type_name["Attractor"] = new FieldFactoryT<AttractorField>();
      map_type_name["AttractorAnisoCurve"] = new FieldFactoryT<AttractorAnisoCurveField>();
    #endif
      map_type_name["MaxEigenHessian"] = new FieldFactoryT<MaxEigenHessianField>();
      background_field = -1;
      boundaryLayer_field = -1;
    }
    
    FieldManager::~FieldManager()
    {
      for(std::map<std::string, FieldFactory*>::iterator it = map_type_name.begin(); 
          it != map_type_name.end(); it++)
        delete it->second;
    }
    
    #if defined(HAVE_POST)
    void Field::putOnNewView()
    {
      if(GModel::current()->getMeshStatus() < 1){
        Msg::Error("No mesh available to create the view: please mesh your model!");
        return;
      }
      std::map<int, std::vector<double> > d;
      std::vector<GEntity*> entities;
      GModel::current()->getEntities(entities);
      for(unsigned int i = 0; i < entities.size(); i++){
        for(unsigned int j = 0; j < entities[i]->mesh_vertices.size(); j++){
          MVertex *v = entities[i]->mesh_vertices[j];
          d[v->getNum()].push_back((*this)(v->x(), v->y(), v->z(), entities[i]));
        }
      }
      std::ostringstream oss;
      oss << "Field " << id;
      PView *view = new PView(oss.str(), "NodeData", GModel::current(), d);
      view->setChanged(true);
    }
    
    void Field::putOnView(PView *view, int comp)
    {
      PViewData *data = view->getData();
      for(int ent = 0; ent < data->getNumEntities(0); ent++){
        for(int ele = 0; ele < data->getNumElements(0, ent); ele++){
          if(data->skipElement(0, ent, ele)) continue;
          for(int nod = 0; nod < data->getNumNodes(0, ent, ele); nod++){
            double x, y, z;
            data->getNode(0, ent, ele, nod, x, y, z);
            double val = (*this)(x, y, z);
            for(int comp = 0; comp < data->getNumComponents(0, ent, ele); comp++)
              data->setValue(0, ent, ele, nod, comp, val);
          }
        }
      }
      std::ostringstream oss;
      oss << "Field " << id;
      data->setName(oss.str());
      data->finalize();
      view->setChanged(true);
      data->destroyAdaptiveData();
    }
    #endif
    
    void FieldManager::setBackgroundMesh(int iView)
    {
      int id = newId();
      Field *f = newField(id, "PostView");
      f->options["IView"]->numericalValue(iView);
      (*this)[id] = f;
      background_field = id;
    }