diff --git a/Plugin/MathEval.cpp b/Plugin/MathEval.cpp
index 6f644249307e510d7fe786140e14a621a24daf5a..b3ba400e09a521f637a89b2e938447b5d5222da8 100644
--- a/Plugin/MathEval.cpp
+++ b/Plugin/MathEval.cpp
@@ -8,6 +8,7 @@
 #include "MathEval.h"
 #include "mathEvaluator.h"
 #include "OctreePost.h"
+#include "GEntity.h"
 
 StringXNumber MathEvalOptions_Number[] = {
   {GMSH_FULLRC, "TimeStep", NULL, -1.},
@@ -96,10 +97,10 @@ PView *GMSH_MathEvalPlugin::execute(PView *view)
   int otherTimeStep = (int)MathEvalOptions_Number[2].def;
   int iOtherView = (int)MathEvalOptions_Number[3].def;
   int forceInterpolation = (int)MathEvalOptions_Number[4].def;
-  int region = (int)MathEvalOptions_Number[5].def;
+  int physicalRegion = (int)MathEvalOptions_Number[5].def;
   std::vector<std::string> expr(9);
   for(int i = 0; i < 9; i++) expr[i] = MathEvalOptions_String[i].def;
-  
+
   PView *v1 = getView(iView, view);
   if(!v1) return view;
   PViewData *data1 = getPossiblyAdaptiveData(v1);
@@ -143,7 +144,7 @@ PView *GMSH_MathEvalPlugin::execute(PView *view)
   }
 
   int numComp2;
-  if(expr[3].size() || expr[4].size() || expr[5].size() || 
+  if(expr[3].size() || expr[4].size() || expr[5].size() ||
      expr[6].size() || expr[7].size() || expr[8].size()){
     numComp2 = 9;
     for(int i = 0; i < 9; i++)
@@ -159,7 +160,7 @@ PView *GMSH_MathEvalPlugin::execute(PView *view)
   }
   expr.resize(numComp2);
 
-  const char *names[] = 
+  const char *names[] =
     { "x", "y", "z", "v0", "v1", "v2", "v3", "v4", "v5", "v6", "v7", "v8",
       "w0", "w1", "w2", "w3", "w4", "w5", "w6", "w7", "w8", "w9" };
   unsigned int numVariables = sizeof(names) / sizeof(names[0]);
@@ -168,7 +169,7 @@ PView *GMSH_MathEvalPlugin::execute(PView *view)
   mathEvaluator f(expr, variables);
   if(expr.empty()) return view;
   std::vector<double> values(numVariables), res(numComp2);
-          
+
   PView *v2 = new PView();
   PViewDataList *data2 = getDataList(v2);
 
@@ -185,7 +186,16 @@ PView *GMSH_MathEvalPlugin::execute(PView *view)
   int timeBeg = (timeStep < 0) ? firstNonEmptyStep : timeStep;
   int timeEnd = (timeStep < 0) ? -timeStep : timeStep + 1;
   for(int ent = 0; ent < data1->getNumEntities(timeBeg); ent++){
-    if (region>0 && (ent+1)!=region) continue;
+    bool ok = (physicalRegion <= 0);
+    if(physicalRegion > 0){
+      GEntity *ge = data1->getEntity(timeBeg, ent);
+      if(ge){
+        std::vector<int>::iterator it = std::find
+          (ge->physicals.begin(), ge->physicals.end(), physicalRegion);
+        ok = (it != ge->physicals.end());
+      }
+    }
+    if(!ok) continue;
     for(int ele = 0; ele < data1->getNumElements(timeBeg, ent); ele++){
       if(data1->skipElement(timeBeg, ent, ele)) continue;
       int numNodes = data1->getNumNodes(timeBeg, ent, ele);
@@ -198,9 +208,9 @@ PView *GMSH_MathEvalPlugin::execute(PView *view)
       std::vector<double> x(numNodes), y(numNodes), z(numNodes);
       for(int nod = 0; nod < numNodes; nod++)
         data1->getNode(timeBeg, ent, ele, nod, x[nod], y[nod], z[nod]);
-      for(int nod = 0; nod < numNodes; nod++) out->push_back(x[nod]); 
-      for(int nod = 0; nod < numNodes; nod++) out->push_back(y[nod]); 
-      for(int nod = 0; nod < numNodes; nod++) out->push_back(z[nod]); 
+      for(int nod = 0; nod < numNodes; nod++) out->push_back(x[nod]);
+      for(int nod = 0; nod < numNodes; nod++) out->push_back(y[nod]);
+      for(int nod = 0; nod < numNodes; nod++) out->push_back(z[nod]);
       for(int step = timeBeg; step < timeEnd; step++){
 	if(!data1->hasTimeStep(step)) continue;
         int step2 = (otherTimeStep < 0) ? step : otherTimeStep;
@@ -228,7 +238,7 @@ PView *GMSH_MathEvalPlugin::execute(PView *view)
       }
     }
   }
-  
+
   if(timeStep < 0){
     for(int i = firstNonEmptyStep; i < data1->getNumTimeSteps(); i++) {
       if(!data1->hasTimeStep(i)) continue;
@@ -237,7 +247,7 @@ PView *GMSH_MathEvalPlugin::execute(PView *view)
   }
   else
     data2->Time.push_back(data1->getTime(timeStep));
-        
+
   data2->setName(data1->getName() + "_MathEval");
   data2->setFileName(data1->getName() + "_MathEval.pos");
   data2->finalize();
diff --git a/Post/PViewData.cpp b/Post/PViewData.cpp
index 2ead5dc066e5b0dbac32037fadb4981f6f0a3b8f..d87e86186c10f77005decd71e61078c77b26fb6b 100644
--- a/Post/PViewData.cpp
+++ b/Post/PViewData.cpp
@@ -27,7 +27,7 @@ PViewData::~PViewData()
 }
 
 bool PViewData::finalize(bool computeMinMax, const std::string &interpolationScheme)
-{ 
+{
   _dirty = false;
   return true;
 }
@@ -95,13 +95,25 @@ void PViewData::setValue(int step, int ent, int ele, int nod, int comp, double v
   Msg::Error("Cannot change field value in this view");
 }
 
+GModel *PViewData::getModel(int step)
+{
+  Msg::Error("Cannot get entity from this view");
+  return 0;
+}
+
+GEntity *PViewData::getEntity(int step, int ent)
+{
+  Msg::Error("Cannot get entity from this view");
+  return 0;
+}
+
 MElement *PViewData::getElement(int step, int ent, int ele)
 {
   Msg::Error("Cannot get element from this view");
-  return NULL;
+  return 0;
 }
 
-void PViewData::setInterpolationMatrices(int type, 
+void PViewData::setInterpolationMatrices(int type,
                                          const fullMatrix<double> &coefVal,
                                          const fullMatrix<double> &expVal)
 {
@@ -110,9 +122,9 @@ void PViewData::setInterpolationMatrices(int type,
   _interpolation[type].push_back(new fullMatrix<double>(expVal));
 }
 
-void PViewData::setInterpolationMatrices(int type, 
+void PViewData::setInterpolationMatrices(int type,
                                          const fullMatrix<double> &coefVal,
-                                         const fullMatrix<double> &expVal, 
+                                         const fullMatrix<double> &expVal,
                                          const fullMatrix<double> &coefGeo,
                                          const fullMatrix<double> &expGeo)
 {
@@ -133,8 +145,8 @@ int PViewData::getInterpolationMatrices(int type, std::vector<fullMatrix<double>
 }
 
 bool PViewData::haveInterpolationMatrices(int type)
-{ 
-  if(!type) 
+{
+  if(!type)
     return !_interpolation.empty();
   else
     return _interpolation.count(type) ? true : false;
@@ -164,39 +176,39 @@ void PViewData::smooth()
 }
 
 bool PViewData::combineTime(nameData &nd)
-{ 
+{
   Msg::Error("Combine time is not implemented for this type of data");
-  return false; 
+  return false;
 }
 
 bool PViewData::combineSpace(nameData &nd)
-{ 
+{
   Msg::Error("Combine space is not implemented for this type of data");
-  return false; 
+  return false;
 }
 
-bool PViewData::searchScalar(double x, double y, double z, double *values, 
+bool PViewData::searchScalar(double x, double y, double z, double *values,
                              int step, double *size)
 {
   if(!_octree) _octree = new OctreePost(this);
   return _octree->searchScalar(x, y, z, values, step, size);
 }
 
-bool PViewData::searchScalarWithTol(double x, double y, double z, double *values, 
+bool PViewData::searchScalarWithTol(double x, double y, double z, double *values,
                                     int step, double *size, double tol)
 {
   if(!_octree) _octree = new OctreePost(this);
   return _octree->searchScalarWithTol(x, y, z, values, step, size, tol);
 }
 
-bool PViewData::searchVector(double x, double y, double z, double *values, 
+bool PViewData::searchVector(double x, double y, double z, double *values,
                              int step, double *size)
 {
   if(!_octree) _octree = new OctreePost(this);
   return _octree->searchVector(x, y, z, values, step, size);
 }
 
-bool PViewData::searchTensor(double x, double y, double z, double *values, 
+bool PViewData::searchTensor(double x, double y, double z, double *values,
                              int step, double *size)
 {
   if(!_octree) _octree = new OctreePost(this);
diff --git a/Post/PViewData.h b/Post/PViewData.h
index e8037e8b592780d49d58c0e267d341e779449c8a..b1130566f0db532b81abf134cd43e18b00aaa54e 100644
--- a/Post/PViewData.h
+++ b/Post/PViewData.h
@@ -77,7 +77,7 @@ class PViewData {
   virtual double getTime(int step){ return 0.; }
 
   // get/set min/max for given step (global over all steps if step=-1)
-  virtual double getMin(int step=-1, bool onlyVisible=false, 
+  virtual double getMin(int step=-1, bool onlyVisible=false,
                         int forceNumComponents=0, int componentMap[9]=0) = 0;
   virtual double getMax(int step=-1, bool onlyVisible=false,
                         int forceNumComponents=0, int componentMap[9]=0) = 0;
@@ -121,7 +121,7 @@ class PViewData {
   // get/set the coordinates and tag of the nod-th node from the
   // ele-th element in the ent-th entity (if the node has a tag,
   // getNode returns it)
-  virtual int getNode(int step, int ent, int ele, int nod, 
+  virtual int getNode(int step, int ent, int ele, int nod,
                       double &x, double &y, double &z){ return 0; }
   virtual void setNode(int step, int ent, int ele, int nod,
                        double x, double y, double z);
@@ -162,9 +162,9 @@ class PViewData {
   virtual int getNumStrings3D(){ return 0; }
 
   // return the i-th 2D/3D string in the view
-  virtual void getString2D(int i, int step, std::string &str, 
+  virtual void getString2D(int i, int step, std::string &str,
                            double &x, double &y, double &style){}
-  virtual void getString3D(int i, int step, std::string &str, 
+  virtual void getString3D(int i, int step, std::string &str,
                            double &x, double &y, double &z, double &style){}
 
   // change the orientation of the ele-th element
@@ -195,13 +195,13 @@ class PViewData {
   adaptiveData *getAdaptiveData(){ return _adaptive; }
 
   // set/get the interpolation matrices for elements with type "type"
-  void setInterpolationMatrices(int type, 
+  void setInterpolationMatrices(int type,
                                 const fullMatrix<double> &coefVal,
                                 const fullMatrix<double> &expVal);
-  void setInterpolationMatrices(int type, 
+  void setInterpolationMatrices(int type,
                                 const fullMatrix<double> &coefVal,
                                 const fullMatrix<double> &expVal,
-                                const fullMatrix<double> &coefGeo, 
+                                const fullMatrix<double> &coefGeo,
                                 const fullMatrix<double> &expGeo);
   int getInterpolationMatrices(int type, std::vector<fullMatrix<double>*> &p);
   bool haveInterpolationMatrices(int type=0);
@@ -217,11 +217,17 @@ class PViewData {
   // combine time steps or elements from multiple datasets
   virtual bool combineTime(nameData &nd);
   virtual bool combineSpace(nameData &nd);
-  
+
   // ask to fill vertex arrays remotely
   virtual bool isRemote(){ return false; }
   virtual int fillRemoteVertexArrays(std::string &options){ return 0; }
 
+  // get GModel (if view supports it)
+  virtual GModel *getModel(int step);
+
+  // get GEntity (if view supports it)
+  virtual GEntity *getEntity(int step, int entity);
+
   // get MElement (if view supports it)
   virtual MElement *getElement(int step, int entity, int element);
 
@@ -230,13 +236,13 @@ class PViewData {
   // post element. If several time steps are present, they are all
   // interpolated unless time step is set to a different value than
   // -1.
-  bool searchScalar(double x, double y, double z, double *values, 
+  bool searchScalar(double x, double y, double z, double *values,
                     int step=-1, double *size=0);
-  bool searchScalarWithTol(double x, double y, double z, double *values, 
+  bool searchScalarWithTol(double x, double y, double z, double *values,
                            int step=-1, double *size=0, double tol=1.e-2);
-  bool searchVector(double x, double y, double z, double *values, 
+  bool searchVector(double x, double y, double z, double *values,
                     int step=-1, double *size=0);
-  bool searchTensor(double x, double y, double z, double *values, 
+  bool searchTensor(double x, double y, double z, double *values,
                     int step=-1, double *size=0);
 
   // I/O routines
diff --git a/Post/PViewDataGModel.cpp b/Post/PViewDataGModel.cpp
index 99ec19ff606ada1152b3f7098ba7ef0943040a9d..59b0e56a54a3517a6cb7b4dee8a8ab4d79616192 100644
--- a/Post/PViewDataGModel.cpp
+++ b/Post/PViewDataGModel.cpp
@@ -34,7 +34,7 @@ static MElement *_getOneElementOfGivenType(GModel *m, int type)
       if((*it)->points.size()) return (*it)->points[0];
     }
     break;
-  case TYPE_LIN: 
+  case TYPE_LIN:
     for(GModel::eiter it = m->firstEdge(); it != m->lastEdge(); it++){
       if((*it)->lines.size()) return (*it)->lines[0];
     }
@@ -126,15 +126,15 @@ bool PViewDataGModel::finalize(bool computeMinMax, const std::string &interpolat
   if(!haveInterpolationMatrices()){
 
     GModel *model = _steps[0]->getModel();
-    
+
     // if an interpolation scheme is explicitly provided, use it
     if(interpolationScheme.size()){
       interpolationMatrices m = _interpolationSchemes[interpolationScheme];
       if(m.size())
-        Msg::Info("Setting interpolation matrices from scheme '%s'", 
+        Msg::Info("Setting interpolation matrices from scheme '%s'",
                   interpolationScheme.c_str());
       else
-        Msg::Error("Could not find interpolation scheme '%s'", 
+        Msg::Error("Could not find interpolation scheme '%s'",
                    interpolationScheme.c_str());
       for(interpolationMatrices::iterator it = m.begin(); it != m.end(); it++){
         if(it->second.size() == 2){
@@ -163,7 +163,7 @@ bool PViewDataGModel::finalize(bool computeMinMax, const std::string &interpolat
                      (int)it->second.size(), interpolationScheme.c_str());
       }
     }
-    
+
     // if we don't have interpolation matrices for a given element
     // type, assume isoparametric elements (except for ElementData,
     // for which we know the interpolation: it's constant)
@@ -199,7 +199,7 @@ bool PViewDataGModel::finalize(bool computeMinMax, const std::string &interpolat
     }
 
   }
-  
+
   return PViewData::finalize();
 }
 
@@ -249,7 +249,7 @@ double PViewDataGModel::getMin(int step, bool onlyVisible, int forceNumComponent
         if(skipElement(step, ent, ele, onlyVisible)) continue;
         for(int nod = 0; nod < getNumNodes(step, ent, ele); nod++){
           double val;
-          getScalarValue(step, ent, ele, nod, val, 
+          getScalarValue(step, ent, ele, nod, val,
                          forceNumComponents, componentMap);
           vmin = std::min(vmin, val);
         }
@@ -275,7 +275,7 @@ double PViewDataGModel::getMax(int step, bool onlyVisible, int forceNumComponent
         if(skipElement(step, ent, ele, onlyVisible)) continue;
         for(int nod = 0; nod < getNumNodes(step, ent, ele); nod++){
           double val;
-          getScalarValue(step, ent, ele, nod, val, 
+          getScalarValue(step, ent, ele, nod, val,
                          forceNumComponents, componentMap);
           vmax = std::max(vmax, val);
         }
@@ -441,6 +441,11 @@ int PViewDataGModel::getNumElements(int step, int ent)
   return _steps[step]->getEntity(ent)->getNumMeshElements();
 }
 
+GEntity *PViewDataGModel::getEntity(int step, int ent)
+{
+  return _steps[step]->getEntity(ent);
+}
+
 MElement *PViewDataGModel::getElement(int step, int ent, int element)
 {
   if(_steps.empty()) return 0;
@@ -638,7 +643,7 @@ void PViewDataGModel::smooth()
     GModel *m = _steps[step]->getModel();
     int numComp = _steps[step]->getNumComponents();
     _steps2.push_back(new stepData<double>(m, numComp, _steps[step]->getFileName(),
-                                           _steps[step]->getFileIndex(), 
+                                           _steps[step]->getFileIndex(),
                                            _steps[step]->getTime()));
     std::map<int, int> nodeConnect;
     for(int ent = 0; ent < getNumEntities(step); ent++){
@@ -770,11 +775,6 @@ bool PViewDataGModel::hasModel(GModel *model, int step)
   return (model == _steps[step]->getModel());
 }
 
-GEntity *PViewDataGModel::getEntity(int step, int ent)
-{
-  return _steps[step]->getEntity(ent);
-}
-
 bool PViewDataGModel::getValueByIndex(int step, int dataIndex, int nod, int comp, double &val)
 {
   double *d = _steps[step]->getData(dataIndex);
diff --git a/Post/PViewDataGModel.h b/Post/PViewDataGModel.h
index 63408da5468948b6fb0309af1df3d133f62bd585..c0a8391c0d4dcc3c95468f781a9ffd1b837fd556 100644
--- a/Post/PViewDataGModel.h
+++ b/Post/PViewDataGModel.h
@@ -46,9 +46,9 @@ class stepData{
   // a set of all "partitions" encountered in the data
   std::set<int> _partitions;
  public:
-  stepData(GModel *model, int numComp, std::string fileName="", int fileIndex=-1, 
+  stepData(GModel *model, int numComp, std::string fileName="", int fileIndex=-1,
            double time=0., double min=VAL_INF, double max=-VAL_INF)
-    : _model(model), _fileName(fileName), _fileIndex(fileIndex), _time(time), 
+    : _model(model), _fileName(fileName), _fileIndex(fileIndex), _time(time),
       _min(min), _max(max), _numComp(numComp), _data(0)
   {
     _model->getEntities(_entities);
@@ -108,7 +108,7 @@ class stepData{
     return _data->size();
   }
   void resizeData(int n)
-  {  
+  {
     if(!_data) _data = new std::vector<real*>(n, (real*)0);
     if(n > (int)_data->size()) _data->resize(n, (real*)0);
   }
@@ -220,17 +220,14 @@ class PViewDataGModel : public PViewData {
   bool hasMultipleMeshes();
   bool hasModel(GModel *model, int step=-1);
   bool useGaussPoints(){ return _type == GaussPointData; }
+  GModel* getModel(int step){ return _steps[step]->getModel(); }
+  GEntity *getEntity(int step, int ent);
+  MElement *getElement(int step, int entity, int element);
 
   // get the data type
   DataType getType(){ return _type; }
-  // direct access to GModel entities
-  GEntity *getEntity(int step, int ent);
   // direct access to value by index
   bool getValueByIndex(int step, int dataIndex, int node, int comp, double &val);
-  // get underlying model
-  GModel* getModel(int step){ return _steps[step]->getModel(); }
-  // get MElement
-  MElement *getElement(int step, int entity, int element);
 
   // Add some data "on the fly" (data is stored in a map, indexed by
   // node or element number depending on the type of dataset)
@@ -238,8 +235,8 @@ class PViewDataGModel : public PViewData {
                int step, double time, int partition, int numComp);
 
   // I/O routines
-  bool readMSH(std::string fileName, int fileIndex, FILE *fp, bool binary, 
-               bool swap, int step, double time, int partition, 
+  bool readMSH(std::string fileName, int fileIndex, FILE *fp, bool binary,
+               bool swap, int step, double time, int partition,
                int numComp, int numNodes, const std::string &interpolationScheme);
   bool writeMSH(std::string fileName, bool binary=false, bool savemesh=true);
   bool readMED(std::string fileName, int fileIndex);