diff --git a/Mesh/Field.cpp b/Mesh/Field.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..8f257c75ad6c2a5131514c4ba64a667acaa338ac
--- /dev/null
+++ b/Mesh/Field.cpp
@@ -0,0 +1,330 @@
+#include<fstream>
+#include "math.h"
+#include "Field.h"
+#include "Gmsh.h"
+#include "Context.h"
+#include "GeoInterpolation.h"
+#include "matheval.h"
+#ifdef HAVE_MATH_EVAL
+#include "matheval.h"
+#endif
+#define MAX_LC 1e22
+extern Context_T CTX;
+FieldManager fields;
+
+void FieldManager::reset(){
+ for(std::map<int,Field*>::iterator it=id_map.begin();it!=id_map.end();it++){
+ delete it->second;
+ }
+ id_map.clear();
+}
+
+Field *FieldManager::get(int id){
+ std::map<int,Field*>::iterator it=id_map.find(id);
+ if(it==id_map.end()){
+ Msg(GERROR,"Field id %i does not exist.\n",id);
+ return NULL;
+ }
+ return it->second;
+}
+int FieldManager::insert(Field *field,int id){
+ if(id==-1){ // get an automatic negative id starting from -1
+ if(id_map.begin()!=id_map.end()){
+ id=id_map.begin()->first-1;
+ }
+ }
+ if(id>0){
+ if(id_map.find(id)!=id_map.end()){
+ Msg(GERROR, "Field id %i is already defined, it will be deleted.\n");
+ delete id_map[id];
+ }
+ id_map[id]=field;
+ return id;
+ }
+}
+
+/* StructuredField*/
+StructuredField::StructuredField(const char *filename){
+ std::ifstream input(filename);
+ 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];
+ data=new double[nt];
+ input.read((char*)data,nt*sizeof(double));
+ input.close();
+}
+StructuredField::~StructuredField(){
+ delete []data;
+}
+double StructuredField::operator()(double x,double y,double z){
+ //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]);
+ 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[j][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;
+}
+
+/*LatLonField*/
+double LatLonField::operator()(double x,double y,double z){
+ return (*field)(asin(z/sqrt(x*x+y*y+z*z)),atan2(y,x),0);
+}
+
+/*ThresholdField*/
+ThresholdField::ThresholdField(Field *_field,
+ double _dmin,double _dmax,double _lcmin,double _lcmax):
+ field(_field),dmin(_dmin),dmax(_dmax),lcmin(_lcmin),lcmax(_lcmax){
+}
+double ThresholdField::operator()(double x,double y,double z){
+ double r=((*field)(x,y,z)-dmin)/(dmax-dmin);
+ r=std::max(std::min(r,1.),0.);
+ double lc=lcmin*(1-r)+lcmax*r;
+ return lc;
+}
+
+/*GradField*/
+GradField::GradField(Field *_field,int _kind,double _delta):
+ field(_field),kind(_kind),delta(_delta){
+ if(delta<0) delta=CTX.lc/10000;
+}
+double GradField::operator()(double x,double y,double z){
+ 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 : /* max */
+ 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(GERROR, "Unknown kind (%i) for GradField.",kind);
+ return 0;
+ }
+}
+
+/*ParametricField*/
+ParametricField::ParametricField(Field *_field,const char *strX,const char *strY,const char *strZ):field(_field){
+ char *sx=strdup(strX),*sy=strdup(strY),*sz=strdup(strZ);
+#if !defined(HAVE_MATH_EVAL)
+ Msg(GERROR, "MathEval is not compiled in this version of Gmsh");
+#else
+ evalX = evaluator_create(sx);
+ evalY = evaluator_create(sy);
+ evalZ = evaluator_create(sz);
+#endif
+}
+ParametricField::~ParametricField(){
+#if !defined(HAVE_MATH_EVAL)
+ Msg(GERROR, "MathEval is not compiled in this version of Gmsh");
+#else
+ evaluator_destroy(evalX);
+ evaluator_destroy(evalY);
+ evaluator_destroy(evalZ);
+#endif
+}
+double ParametricField::operator()(double x,double y,double z){
+#if !defined(HAVE_MATH_EVAL)
+ Msg(GERROR, "MathEval is not compiled in this version of Gmsh");
+#else
+ static char *names[3] = {"x","y","z"};
+ double values [3] = {x,y,z};
+ const double nx = evaluator_evaluate(evalX, 3, names, values);
+ const double ny = evaluator_evaluate(evalY, 3, names, values);
+ const double nz = evaluator_evaluate(evalZ, 3, names, values);
+ return (*field)(x,y,z);
+#endif
+}
+
+/* FunctionField */
+FunctionField::FunctionField(std::list<Field*> *_list,const char *str):list(_list){
+#if !defined(HAVE_MATH_EVAL)
+ Msg(GERROR, "MathEval is not compiled in this version of Gmsh");
+#else
+ char *s=strdup(str);
+ eval = evaluator_create(s);
+ values=new double[3+list->size()];
+ names=new char*[3+list->size()];
+ names[0]=strdup("x");
+ names[1]=strdup("y");
+ names[2]=strdup("z");
+ int p=3;
+ for(std::list<Field*>::iterator it=list->begin();it!=list->end();it++){
+ asprintf(names+p,"f%i",p-3);
+ p++;
+ }
+#endif
+}
+double FunctionField::operator()(double x,double y,double z){
+#if !defined(HAVE_MATH_EVAL)
+ Msg(GERROR, "MathEval is not compiled in this version of Gmsh");
+#else
+ values[0]=x;
+ values[1]=y;
+ values[2]=z;
+ int p=3;
+ for(std::list<Field*>::iterator it=list->begin();it!=list->end();it++){
+ values[p]=(**it)(x,y,z);
+ p++;
+ }
+ return evaluator_evaluate(eval,p,names,values);
+#endif
+}
+FunctionField::~FunctionField(){
+#if !defined(HAVE_MATH_EVAL)
+ Msg(GERROR, "MathEval is not compiled in this version of Gmsh");
+#else
+ int n=3+list->size();
+ for(int i=0;i<n;i++){
+ free(names[i]);
+ }
+ delete []names;
+ delete []values;
+ delete list;
+ evaluator_destroy(eval);
+#endif
+}
+
+/* PostViewField */
+double PostViewField::operator()(double x, double y, double z) {
+ double l = 0.;
+ double fact[9] = {0.001, 0.0025, 0.005, 0.0075, 0.01, 0.025, 0.05, 0.075, 0.1};
+ if(!octree->searchScalar(x, y, z, &l, 0)){
+ // try really hard to find an element around the point
+ for(int i = 0; i < 9; i++){
+ double eps = CTX.lc * fact[i];
+ if(octree->searchScalar(x + eps, y, z, &l, 0)) break;
+ if(octree->searchScalar(x - eps, y, z, &l, 0)) break;
+ if(octree->searchScalar(x, y + eps, z, &l, 0)) break;
+ if(octree->searchScalar(x, y - eps, z, &l, 0)) break;
+ if(octree->searchScalar(x, y, z + eps, &l, 0)) break;
+ if(octree->searchScalar(x, y, z - eps, &l, 0)) break;
+ if(octree->searchScalar(x + eps, y - eps, z - eps, &l, 0)) break;
+ if(octree->searchScalar(x + eps, y + eps, z - eps, &l, 0)) break;
+ if(octree->searchScalar(x - eps, y - eps, z - eps, &l, 0)) break;
+ if(octree->searchScalar(x - eps, y + eps, z - eps, &l, 0)) break;
+ if(octree->searchScalar(x + eps, y - eps, z + eps, &l, 0)) break;
+ if(octree->searchScalar(x + eps, y + eps, z + eps, &l, 0)) break;
+ if(octree->searchScalar(x - eps, y - eps, z + eps, &l, 0)) break;
+ if(octree->searchScalar(x - eps, y + eps, z + eps, &l, 0)) break;
+ }
+ }
+ if(l <= 0) return MAX_LC;
+ return l;
+}
+
+PostViewField::PostViewField(Post_View *view){
+ Msg(INFO, "Field from '%s'.", view->Name);
+ octree = new OctreePost(view);
+ view_num = view->Num;
+}
+
+/* Min Field */
+MinField::MinField(){};
+double MinField::operator()(double x,double y,double z){
+ if(size()==0){
+ Msg(GERROR, "Requesting minimum of a void list of field.");
+ return 0;
+ }
+ iterator it=begin();
+ double v=(**it++)(x,y,z);
+ while(it!=end()){
+ v=std::min(v,(**it++)(x,y,z));
+ }
+ return v;
+}
+
+/* Attractor Field */
+#define maxpts 1
+void AttractorField::addPoint ( double X, double Y, double Z, double lc){
+ attractorPoints.push_back(SPoint3(X,Y,Z));
+}
+AttractorField::~AttractorField(){
+#ifdef HAVE_ANN_
+ if(kdtree) delete kdtree;
+ if(zeronodes) annDeallocPts(zeronodes);
+ delete [] index;
+ delete [] dist;
+#endif
+}
+AttractorField::AttractorField()
+#ifdef HAVE_ANN_
+ : kdtree (0), zeronodes(0)
+#endif
+{
+#ifdef HAVE_ANN_
+ index = new ANNidx[maxpts];
+ dist = new ANNdist[maxpts];
+#endif
+}
+
+void AttractorField::buildFastSearchStructures(){
+#ifdef HAVE_ANN_
+ if(zeronodes){
+ annDeallocPts(zeronodes);
+ delete kdtree;
+ }
+ int totpoints = attractorPoints.size();
+ if (totpoints)
+ zeronodes = annAllocPts(totpoints, 4);
+ int k = 0;
+ for (std::list <SPoint3>::iterator it2 = attractorPoints.begin();
+ it2 != attractorPoints.end(); ++it2){
+ zeronodes[k][0]=it2->x();
+ zeronodes[k][1]=it2->y();
+ zeronodes[k++][2]=it2->z();
+ }
+ kdtree = new ANNkd_tree(zeronodes, totpoints, 3);
+#endif
+}
+
+double AttractorField::operator()(double X, double Y, double Z){
+ double xyz[3] = {X, Y, Z};
+ if(attractorPoints.size() == 1){
+ SPoint3 p = *(attractorPoints.begin());
+ return sqrt((p.x()-X)*(p.x()-X)+(p.y()-Y)*(p.y()-Y)+(p.z()-Z)*(p.z()-Z));
+ }
+ else{
+#ifdef HAVE_ANN_
+ kdtree->annkSearch(xyz, maxpts, index, dist);
+ return sqrt(dist[0]);
+#else
+ Msg(GERROR,"GMSH should be compiled with ANN in order to enable attractors !");
+#endif
+ }
+}
+void AttractorField::addCurve(Curve *c, int N) {
+ for(int i = 0; i < N; i++){
+ double u = (double)i / (N - 1);
+ Vertex V = InterpolateCurve (c, u, 0);
+ addPoint(V.Pos.X, V.Pos.Y, V.Pos.Z);
+ }
+}
+void AttractorField::addGEdge(GEdge *c, int N) {
+ for(int i = 0; i < N; i++){
+ double u = (double)i / (N - 1);
+ Range<double> b = c->parBounds(0);
+ double t = b.low() + u * (b.high() - b.low());
+ GPoint gp = c->point(t);
+ addPoint(gp.x(), gp.y(), gp.z());
+ }
+}
diff --git a/Mesh/Field.h b/Mesh/Field.h
new file mode 100644
index 0000000000000000000000000000000000000000..a4252db6aa2bcf61eb50f3eaa6d89a80b3bf85cf
--- /dev/null
+++ b/Mesh/Field.h
@@ -0,0 +1,105 @@
+#ifndef FIELD_H_
+#define FIELD_H_
+#include<map>
+#include<list>
+#include "Views.h"
+#ifdef HAVE_ANN_
+#include "ANN/ANN.h"
+#endif
+#include "Geo.h"
+#include "GEdge.h"
+#include "OctreePost.h"
+class Field{
+ public:
+ virtual double operator()(double x,double y,double z)=0;
+};
+class FieldManager{
+ std::map<int,Field*> id_map;
+ public:
+ void reset();
+ int insert(Field *field,int id=-1 /* -1= automatic id */); /* return the id */
+ Field *get(int id);
+};
+extern FieldManager fields;
+
+class StructuredField:public Field{
+ double o[3],d[3];
+ int n[3];
+ double *data;
+ public :
+ StructuredField(const char *filename);
+ virtual ~StructuredField();
+ double operator()(double x,double y,double z);
+};
+class LatLonField:public Field{
+ Field *field;
+ public:
+ LatLonField(Field *_field):field(_field){};
+ double operator()(double x,double y,double z);
+};
+class ParametricField:public Field{
+ void *evalX,*evalY,*evalZ;
+ Field *field;
+ public:
+ ParametricField(Field *_field,const char *strX,const char *strY,const char *strZ);
+ double operator()(double x,double y,double z);
+ ~ParametricField();
+};
+class ThresholdField:public Field{
+ Field *field;
+ double dmin,dmax,lcmin,lcmax;
+ public:
+ ThresholdField(Field *field,double dmin,double dmax,double lcmin,double lcmax);
+ double operator()(double x,double y,double z);
+};
+class GradField:public Field{
+ Field *field;
+ int kind;/*0 x, 1 y, 2 z, 3 max */
+ double delta;
+ public:
+ GradField(Field *field,int _kind=2,double _delta=-1);
+ double operator()(double x,double y,double z);
+};
+class FunctionField:public Field{
+ std::list<Field*> *list;
+ char **names;
+ double *values;
+ void *eval;
+ public:
+ FunctionField(std::list<Field*> *_list,const char *fct);
+ double operator()(double x,double y,double z);
+ ~FunctionField();
+};
+class PostViewField:public Field{
+ OctreePost *octree;
+ int view_num;
+ public:
+ PostViewField(Post_View *view);
+ double operator()(double x,double y,double z);
+};
+class MinField:public Field,public std::list<Field*>{
+ public:
+ MinField();
+ double operator()(double x,double y,double z);
+};
+
+class AttractorField:public Field{
+ protected :
+#ifdef HAVE_ANN_
+ ANNkd_tree* kdtree;
+ ANNpointArray zeronodes;
+ ANNidxArray index;
+ ANNdistArray dist;
+#endif
+ std::list <SPoint3> attractorPoints;
+public :
+ AttractorField();
+ virtual ~AttractorField();
+ void addPoint ( double X, double Y, double Z, double lc = -1.0);
+ void buildFastSearchStructures() ;
+ void addGEdge(GEdge*, int);
+ void addCurve(Curve*, int);
+ double operator()(double X, double Y, double Z) ;
+};
+
+#endif