Field.cpp

  {
    return "MathEval";
  }
  std::string getDescription()
  {
    return "Evaluate a mathematical expression. The expression can contain\n"
      "x, y, z for spatial coordinates, F0, F1, ... for field values, and\n"
      "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\n"
      "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) 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_NO_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)
  {
    // FIXME: should test unique view num instead, but that would be slower
    if(view_index < 0 || view_index >= (int)PView::list.size())
      return MAX_LC;
    if(update_needed){
      if(octree) delete octree;
      octree = new OctreePost(PView::list[view_index]);
      update_needed = false;
    }
    double l = 0.;
    if(!octree->searchScalarWithTol(x, y, z, &l, 0, 0, 10.))
      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\n"
       "option is needed for backward compatibility with the BackgroundMesh option",
       &update_needed);
  }
  ~PostViewField()
  {
    if(octree) delete octree;
  }
};
#endif

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) 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) 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\n"
      "curves, surfaces or volumes.";
  }
  double operator() (double x, double y, double z, GEntity *ge=0)
  {
    Field *f = (GModel::current()->getFields()->get(ifield));
    if(!f) 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)
class AttractorField : public Field
{
  ANNkd_tree *kdtree;
  ANNpointArray zeronodes;
  ANNidxArray index;
  ANNdistArray dist;
  std::list<int> nodes_id, edges_id, faces_id;
  int n_nodes_by_edge;
 public:
  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 geomtric 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 discetized each curve",
       &update_needed);
    options["FacesList"] = new FieldOptionList
      (faces_id, "Indices of surfaces in the geometric model (Warning: might\n"
       "give strange results for complex surfaces)", &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\n"
      "be used to compute the distance from curves, in which case each curve\n"
      "is replaced by NNodesByEdge equidistant nodes and the distance from those\n"
      "nodes is computed.";
  }
  virtual double operator() (double X, double Y, double Z, GEntity *ge=0)
  {
    if(update_needed) {
      if(zeronodes) {
        annDeallocPts(zeronodes);
        delete kdtree;
      }
      int totpoints = nodes_id.size() + n_nodes_by_edge * edges_id.size() + 
        n_nodes_by_edge * n_nodes_by_edge * faces_id.size();
      if(totpoints)
        zeronodes = annAllocPts(totpoints, 4);
      int k = 0;
      for(std::list<int>::iterator it = nodes_id.begin();
          it != nodes_id.end(); ++it) {
        Vertex *v = FindPoint(*it);
        if(v) {
          zeronodes[k][0] = v->Pos.X;
          zeronodes[k][1] = v->Pos.Y;
          zeronodes[k++][2] = v->Pos.Z;
        }
        else {
          GVertex *gv = GModel::current()->getVertexByTag(*it);
          if(gv) {
            zeronodes[k][0] = gv->x();
            zeronodes[k][1] = gv->y();
            zeronodes[k++][2] = gv->z();
          }
        }
      }
      for(std::list<int>::iterator it = edges_id.begin();
          it != edges_id.end(); ++it) {
        Curve *c = FindCurve(*it);
        if(c) {
          for(int i = 0; i < n_nodes_by_edge; 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;
          }
        }
        else {
          GEdge *e = GModel::current()->getEdgeByTag(*it);
          if(e) {
            for(int i = 0; i < n_nodes_by_edge; 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);
              zeronodes[k][0] = gp.x();
              zeronodes[k][1] = gp.y();
              zeronodes[k++][2] = gp.z();
            }
          }
        }
      }
      // 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);
              zeronodes[k][0] = V.Pos.X;
              zeronodes[k][1] = V.Pos.Y;
              zeronodes[k++][2] = V.Pos.Z;
            }
          }
        }
        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);
                zeronodes[k][0] = gp.x();
                zeronodes[k][1] = gp.y();
                zeronodes[k++][2] = gp.z();
              }
            }
          }
        }
      }
      kdtree = new ANNkd_tree(zeronodes, totpoints, 3);
      update_needed = false;
    }
    double xyz[3] = { X, Y, Z };
    kdtree->annkSearch(xyz, 1, index, dist);
    return sqrt(dist[0]);
  }
};

#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>();
  map_type_name["BoundaryLayer"] = new FieldFactoryT<BoundaryLayerField>();
  map_type_name["Box"] = new FieldFactoryT<BoxField>();
  map_type_name["Cylinder"] = new FieldFactoryT<CylinderField>();
  map_type_name["LonLat"] = new FieldFactoryT<LonLatField>();
#if !defined(HAVE_NO_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["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>();
#if defined(HAVE_ANN)
  map_type_name["Attractor"] = new FieldFactoryT<AttractorField>();
#endif
  map_type_name["MaxEigenHessian"] = new FieldFactoryT<MaxEigenHessianField>();
  background_field = -1;
}

#if !defined(HAVE_NO_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;
}