add structured grid function

Solver/dgGroupOfElements.cpp

@@ -337,8 +337,10 @@ dgGroupOfFaces::dgGroupOfFaces (const dgGroupOfElements &elGroup, std::string bo
   _groupLeft(elGroup),_groupRight(elGroup)
 {
   _boundaryTag=boundaryTag;
-  if(boundaryTag=="")
-    throw;
+  if(boundaryTag==""){
+    Msg::Warning ("empty boundary tag, group of boundary faces not created");
+    return;
+  }
   _fsLeft=_groupLeft.getElement(0)->getFunctionSpace(pOrder);
   _closuresLeft = _fsLeft->vertexClosure;
   _fsRight=NULL;
@@ -357,8 +359,10 @@ dgGroupOfFaces::dgGroupOfFaces (const dgGroupOfElements &elGroup, std::string bo
   _groupLeft(elGroup),_groupRight(elGroup)
 {
   _boundaryTag=boundaryTag;
-  if(boundaryTag=="")
-    throw;
+  if(boundaryTag==""){
+    Msg::Warning ("empty boundary tag, group of boundary faces not created");
+    return;
+  }
   _fsLeft=_groupLeft.getElement(0)->getFunctionSpace(pOrder);
   _closuresLeft = _fsLeft->edgeClosure;
   _fsRight=NULL;

Solver/function.cpp

+#include <sstream>
+#include <fstream>
 #include "function.h"
 #include "SPoint3.h"
 #include "MElement.h"
-#include <sstream>
 
 // dataCache members
 dataCache::dataCache(dataCacheMap *cacheMap) : _valid(false) {
@@ -37,8 +38,10 @@ function *function::get(std::string functionName, bool acceptNull)
   if (it==_allFunctions.end()) {
     if (acceptNull)
       return NULL;
-    else
+    else{
+      Msg::Error("unknown function : '%s'\n",functionName.c_str());
       throw;
+    }
   }
   return it->second;
 }
@@ -115,6 +118,83 @@ class functionXYZ : public function {
   }
 };
 
+class functionStructuredGridFile : public function {
+  public:
+  class data;
+  std::string _coordFunction;
+  int n[3];
+  double d[3],o[3];
+  double get(int i,int j, int k){
+    return v[(i*n[1]+j)*n[2]+k];
+  }
+  double *v;
+  class data : public dataCacheDouble{
+    dataCacheDouble &coord;
+    functionStructuredGridFile *_fun;
+    public:
+    data(functionStructuredGridFile *fun, dataCacheMap *m):
+      dataCacheDouble(*m,m->getNbEvaluationPoints(),1),
+      coord(m->get(fun->_coordFunction,this))
+    {
+      _fun=fun;
+    }
+    void _eval(){
+      for(int pt=0;pt<_value.size1();pt++){
+        double xi[3];
+        int id[3];
+        for(int i=0;i<3;i++){
+          id[i]=(int)((coord(pt,i)-_fun->o[i])/_fun->d[i]);
+	int a=id[i];
+          id[i]=std::max(0,std::min(_fun->n[i]-1,id[i]));
+          xi[i]=(coord(pt,i)-_fun->o[i]-id[i]*_fun->d[i])/_fun->d[i];
+          xi[i]=std::min(1.,std::max(0.,xi[i]));
+        }
+        _value(pt,0) =
+           _fun->get(id[0]   ,id[1]   ,id[2]   )*(1-xi[0])*(1-xi[1])*(1-xi[2])
+          +_fun->get(id[0]   ,id[1]   ,id[2]+1 )*(1-xi[0])*(1-xi[1])*(  xi[2])
+          +_fun->get(id[0]   ,id[1]+1 ,id[2]   )*(1-xi[0])*(  xi[1])*(1-xi[2])
+          +_fun->get(id[0]   ,id[1]+1 ,id[2]+1 )*(1-xi[0])*(  xi[1])*(  xi[2])
+          +_fun->get(id[0]+1 ,id[1]   ,id[2]   )*(  xi[0])*(1-xi[1])*(1-xi[2])
+          +_fun->get(id[0]+1 ,id[1]   ,id[2]+1 )*(  xi[0])*(1-xi[1])*(  xi[2])
+          +_fun->get(id[0]+1 ,id[1]+1 ,id[2]   )*(  xi[0])*(  xi[1])*(1-xi[2])
+          +_fun->get(id[0]+1 ,id[1]+1 ,id[2]+1 )*(  xi[0])*(  xi[1])*(  xi[2]);
+      }
+    }
+  };
+  dataCacheDouble *newDataCache(dataCacheMap* m) {
+    return new data(this,m);
+  }
+  functionStructuredGridFile(const std::string filename, const std::string coordFunction){
+    _coordFunction=coordFunction;
+    std::ifstream input(filename.c_str());
+    if(!input)
+      Msg::Error("cannot open file : %s",filename.c_str());
+    if(filename.substr(filename.size()-4,4)!=".bin") {
+      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];
+      v = new double [nt];
+      for (int i=0; i<nt; i++){
+        input>>v[i];
+      }
+    } else {
+      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];
+      v = new double[nt];
+      input.read((char *)v, nt * sizeof(double));
+    }
+    static int c=0;
+    std::ostringstream oss;
+    oss<<"FunctionStructured"<<c++;
+    _name = oss.str();
+    function::add(_name,this);
+  }
+  ~functionStructuredGridFile(){
+    delete []v;
+  }
+};
+
 
 // constant values copied over each line
 class functionConstant::data : public dataCacheDouble {
@@ -180,7 +260,6 @@ public:
 };
 */
 
-
 #include "Bindings.h"
 
 void function::registerDefaultFunctions()
@@ -193,11 +272,20 @@ void function::registerBindings(binding *b){
   methodBinding *mb;
   mb = cb->addMethod("getName",&function::getName);
   mb->setDescription("Return the string which identifies this function.");
+
   cb = b->addClass<functionConstant>("functionConstant");
   cb->setDescription("A constant (scalar or vector) function");
   mb = cb->setConstructor<functionConstant,fullMatrix<double>*>();
   mb->setArgNames("v",NULL);
   mb->setDescription("A new constant function wich values 'v' everywhere. v can be a row-vector.");
   cb->setParentClass<function>();
+
+  cb = b->addClass<functionStructuredGridFile>("functionStructuredGridFile");
+  cb->setParentClass<function>();
+  cb->setDescription("A function to interpolate through data given on a structured grid");
+  mb = cb->setConstructor<functionStructuredGridFile,std::string, std::string>();
+  mb->setArgNames("fileName","coordinateFunction",NULL);
+  mb->setDescription("Tri-linearly interpolate through data in file 'fileName' at coordinate given by 'coordinateFunction'.\nThe file format is :\nx0 y0 z0\ndx dy dz\nnx ny nz\nv(0,0,0) v(0,0,1) v(0 0 2) ..."); 
+
 }