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 {