Field.cpp

// Gmsh - Copyright (C) 1997-2013 C. Geuzaine, J.-F. Remacle
//
// See the LICENSE.txt file for license information. Please report all
// bugs and problems to the public mailing list <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"
#include "CenterlineField.h"
#include "STensor3.h"
#include "meshMetric.h"
#if defined(HAVE_POST)
#include "PView.h"
#include "OctreePost.h"
#include "PViewDataList.h"
#include "MVertex.h"
#endif

#if defined(HAVE_ANN)
#include "ANN/ANN.h"
#endif

Field::~Field()
{
  for(std::map<std::string, FieldOption*>::iterator it = options.begin();
      it != options.end(); ++it)
    delete it->second;
  for(std::map<std::string, FieldCallback*>::iterator it = callbacks.begin();
      it != callbacks.end(); ++it)
    delete it->second;
}

FieldOption *Field::getOption(const std::string optionName)
{
  std::map<std::string, FieldOption*>::iterator it = options.find(optionName);
  if (it == options.end()) {
    Msg::Error("field option :%s does not exist", optionName.c_str());
    return NULL;
  }
  return it->second;
}

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, outside_value_set;
  double outside_value;
  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);
    options["SetOutsideValue"] = new FieldOptionBool(outside_value_set, "True to use the \"OutsideValue\" option. If False, the last values of the grid are used.");
    options["OutsideValue"] = new FieldOptionDouble(outside_value, "Value of the field outside the grid (only used if the \"SetOutsideValue\" option is true).");
    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;
      if (outside_value_set && (id[0][i] < 0 || id[1][i] >= n[i]) && n[i] > 1) 
        return outside_value;
      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 FrustumField : public Field
{
  double x1,y1,z1;
  double x2,y2,z2;
  double r1i,r1o,r2i,r2o;
  double v1i,v1o,v2i,v2o;

 public:
  std::string getDescription()
  {
    return "This field is an extended cylinder with inner (i) and outer (o) radiuses"
      "on both endpoints (1 and 2). Length scale is bilinearly interpolated between"
      "these locations (inner and outer radiuses, endpoints 1 and 2)"
      "The field values for a point P are given by :"
      "  u = P1P.P1P2/||P1P2|| "
      "  r = || P1P - u*P1P2 || "
      "  Ri = (1-u)*R1i + u*R2i "
      "  Ro = (1-u)*R1o + u*R2o "
      "  v = (r-Ri)/(Ro-Ri)"
      "  lc = (1-v)*( (1-u)*v1i + u*v2i ) + v*( (1-u)*v1o + u*v2o )"
      "    where (u,v) in [0,1]x[0,1]";
  }
  FrustumField()
  {
    x1 = y1 = z1 = 0.;
    x2 = y2 = 0.;
    z1 = 1.;
    r1i = r2i = 0.;
    r1o = r2o = 1.;
    v1i = v2i = 0.1;
    v1o = v2o = 1.;

    options["X1"] = new FieldOptionDouble
      (x1, "X coordinate of endpoint 1");
    options["Y1"] = new FieldOptionDouble
      (y1, "Y coordinate of endpoint 1");
    options["Z1"] = new FieldOptionDouble
      (z1, "Z coordinate of endpoint 1");
    options["X2"] = new FieldOptionDouble
      (x2, "X coordinate of endpoint 2");
    options["Y2"] = new FieldOptionDouble
      (y2, "Y coordinate of endpoint 2");
    options["Z2"] = new FieldOptionDouble
      (z2, "Z coordinate of endpoint 2");

    options["R1_inner"] = new FieldOptionDouble
      (r1i, "Inner radius of Frustum at endpoint 1");
    options["R1_outer"] = new FieldOptionDouble
      (r1o, "Outer radius of Frustum at endpoint 1");
    options["R2_inner"] = new FieldOptionDouble
      (r2i, "Inner radius of Frustum at endpoint 2");
    options["R2_outer"] = new FieldOptionDouble
      (r2o, "Outer radius of Frustum at endpoint 2");

    options["V1_inner"] = new FieldOptionDouble
      (v1i, "Element size at point 1, inner radius");
    options["V1_outer"] = new FieldOptionDouble
      (v1o, "Element size at point 1, outer radius");
    options["V2_inner"] = new FieldOptionDouble
      (v2i, "Element size at point 2, inner radius");
    options["V2_outer"] = new FieldOptionDouble
      (v2o, "Element size at point 2, outer radius");

  }
  const char *getName()
  {
    return "Frustum";
  }
  double operator() (double x, double y, double z, GEntity *ge=0)
  {
    double dx = x-x1;
    double dy = y-y1;
    double dz = z-z1;
    double x12 = x2-x1;
    double y12 = y2-y1;
    double z12 = z2-z1;
    double l12 = sqrt( x12*x12 + y12*y12 + z12*z12 );

    double l = (dx*x12 + dy*y12 + dz*z12)/l12 ;
    double r = sqrt( dx*dx+dy*dy+dz*dz - l*l );

    double u = l/l12 ;  // u varies between 0 (P1) and 1 (P2)
    double ri = (1-u)*r1i + u*r2i ;
    double ro = (1-u)*r1o + u*r2o ;
    double v = (r-ri)/(ro-ri) ;  // v varies between 0 (inner) and 1 (outer)

    double lc = MAX_LC ;
    if( u>=0 && u<=1 && v>=0 && v<=1 ){
      lc = (1-v)*( (1-u)*v1i + u*v2i ) + v*( (1-u)*v1o + u*v2o );
    }
    return lc;

  }
};

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++){