From 73469ba10b5604cb40125e48479134ff39d88e8f Mon Sep 17 00:00:00 2001
From: vigmond <ed.vigmond@gmail.com>
Date: Wed, 18 Nov 2020 20:50:05 +0100
Subject: [PATCH] Display element based data

---
 src/CutSurfaces.cc    | 10 ++++++----
 src/CutSurfaces.h     |  2 +-
 src/DataAllInMem.h    |  2 +-
 src/DataOnDemand.h    |  4 ++--
 src/DrawingObjects.cc | 13 ++++++++-----
 src/DrawingObjects.h  | 28 ++++++++++++++++------------
 src/Model.cc          | 42 ++++++++++++++++++++++++++++++++++++++----
 src/Model.h           |  1 +
 src/PolyGon.cc        | 10 ++++++----
 src/Quadrilateral.cc  | 26 +++++++++++++++-----------
 src/Surfaces.cc       | 25 +++++++++++++++----------
 src/Surfaces.h        | 21 ++++++++++++---------
 src/TBmeshWin.cc      | 42 ++++++++++++++++++++++++++----------------
 src/ThreadedData.h    |  2 +-
 src/Triangle.cc       | 18 ++++++++++--------
 15 files changed, 158 insertions(+), 88 deletions(-)

diff --git a/src/CutSurfaces.cc b/src/CutSurfaces.cc
index 79d2cde..b06d8b4 100644
--- a/src/CutSurfaces.cc
+++ b/src/CutSurfaces.cc
@@ -18,20 +18,21 @@ CutSurfaces::CutSurfaces(GLfloat *n)
 
 /** add a Surface element to the surface
  *
- *  \param se the element
- *  \param n  the normal
- *  \param ni list for data interpolation
+ *  \param se   the element
+ *  \param ni   list for data interpolation
+ *  \param sele element being cut
  *
  *  \pre each element pointer points to a single element and not a list
  *
  *  \note points are not shared between elements
  */
 void
-CutSurfaces :: addEle( SurfaceElement *se, Interpolator<DATA_TYPE> *ni )
+CutSurfaces :: addEle( SurfaceElement *se, Interpolator<DATA_TYPE> *ni, int sele )
 {
   _ele.push_back(se);
   _ptarr.push_back(const_cast<GLfloat *>(const_cast<PPoint*>(se->pt())->pt()));
   _interp.push_back(ni);
+  _sele.push_back(sele);
 }
 
 
@@ -59,6 +60,7 @@ CutSurfaces::clear( GLfloat *n )
       delete[] _interp[i];
     }
   _ele.clear();
+  _sele.clear();
   _ptarr.clear();
   _interp.clear();
   _zlist.clear();
diff --git a/src/CutSurfaces.h b/src/CutSurfaces.h
index 831257c..ee7ca17 100644
--- a/src/CutSurfaces.h
+++ b/src/CutSurfaces.h
@@ -20,7 +20,7 @@ class CutSurfaces : public Surfaces
     ~CutSurfaces();
     void get_vert_norms( GLfloat *vn ){}
     void determine_vert_norms( PPoint & ){}
-    void addEle( SurfaceElement *e, Interpolator<DATA_TYPE>* );
+    void addEle( SurfaceElement *e, Interpolator<DATA_TYPE>*, int=-1 );
     DATA_TYPE  interpolate( int e, DATA_TYPE *d, int p )
     {
       return _interp[e][p].interpolate( d );
diff --git a/src/DataAllInMem.h b/src/DataAllInMem.h
index e4c7c1b..131643f 100644
--- a/src/DataAllInMem.h
+++ b/src/DataAllInMem.h
@@ -32,7 +32,7 @@ class DataAllInMem : public DataClass<T>
     virtual T         max();	        // maximum value
     virtual T         min(int);         // minimum value at a time
     virtual T         min();            // minimum value
-    virtual T*        slice(int);       // return pointer to data slice
+    virtual T*        slice(int=-1);       // return pointer to data slice
     virtual void      time_series( int, T* );  // time series for a point
     inline virtual void  increment(int){}
 };
diff --git a/src/DataOnDemand.h b/src/DataOnDemand.h
index 362d956..b5b187d 100644
--- a/src/DataOnDemand.h
+++ b/src/DataOnDemand.h
@@ -25,7 +25,7 @@ class DataOnDemand: public DataClass<T>
     virtual T      max(){assert(0);}                    //!< maximum data value at a time
     virtual T      min(int t){return extremum(t,false);}//!< minimum data value at a time
     virtual T      min(){assert(0);}                    //!< minimum data value at a time
-    virtual T*     slice(int);                          //!< pointer to time slice of data
+    virtual T*     slice(int=-1);                       //!< pointer to time slice of data
     virtual void   time_series( int, T* );              //!< time series for a point
     virtual void   increment(int){}                     //!< time slice increment
 
@@ -223,7 +223,7 @@ T* DataOnDemand<T>::slice( int tm )
 {
   if ( tm>maxtm )
     return NULL;
-  if( tm == last_tm )
+  if( tm == last_tm || tm == -1 )
     return data;
 
   if( _ftype == FTIGB || _ftype == FTDynPt ) {
diff --git a/src/DrawingObjects.cc b/src/DrawingObjects.cc
index aec7d24..ed1c9c5 100644
--- a/src/DrawingObjects.cc
+++ b/src/DrawingObjects.cc
@@ -7,10 +7,10 @@
 
 #ifdef _OPENMP
 #include <omp.h>
-omp_lock_t writelock;
-bool init_MPI_lock=false;
-#endif
+static omp_lock_t writelock;
+#endif 
 
+bool MultiPoint::_init_lock = false;
 
 enum lpint_enum{ BOTH_ON_PLANE, NO_INTERSECTION };
 
@@ -126,13 +126,16 @@ void SurfaceElement::read_normals( int e0, int e1, const char *fnb )
 }
 
 
+/**
+ * @brief initialize OMP lock
+ */
 void
 MultiPoint::init_lock()
 {
 #ifdef _OPENMP
-  if( !init_MPI_lock ) {
+  if( !_init_lock ) {
     omp_init_lock(&writelock);
-    init_MPI_lock=true;
+    _init_lock=true;
   }
 #endif 
 }
diff --git a/src/DrawingObjects.h b/src/DrawingObjects.h
index 0d05956..7c38fb1 100644
--- a/src/DrawingObjects.h
+++ b/src/DrawingObjects.h
@@ -49,6 +49,10 @@ gzFile openFile( const char*, const char* );
 #define MAX_NUM_SURF 6          // per element
 #define MAX_NUM_SURF_NODES 4    // per element
 
+// colouring modes
+#define CM_FLAT 1
+#define CM_ELEM 2
+
 class DrawingObj
 {
   public:
@@ -138,7 +142,8 @@ class MultiPoint : public DrawingObj
     PPoint* _pt;            //!< pointer to point list
     int    _ptsPerObj;		//!< \#nodes to define one object
     int    _nedge;          //!< \#edges
-    void   init_lock();
+    void   init_lock();     //!< initialize OMP lock
+    static bool _init_lock; //!< is lock initialized?
 };
 
 
@@ -177,16 +182,15 @@ class SurfaceElement : public MultiPoint
     SurfaceElement(PPoint *p, int n):MultiPoint(p,n,n),_nrml(NULL),_ptnrml(NULL) {}
     virtual ~SurfaceElement(){if(_nrml) delete[] _nrml;if(_ptnrml)delete[] _ptnrml;}
 
-    virtual void     compute_normals( int, int )=0;
-    const   void     nrml( GLfloat *n ){ _nrml=n; };
+    virtual void  compute_normals( int, int )=0;
+    const   void  nrml( GLfloat *n ){ _nrml=n; };
   inline  const   GLfloat* nrml( int a=0 ) {return _nrml==NULL?NULL:_nrml+3*a; }
   inline  const   GLfloat* ptnrml( int a=0 ) {return _ptnrml==NULL?NULL:_ptnrml+3*a; }
-    virtual int      buffer( int, GLfloat*, Colourscale*, DATA_TYPE*,
-                         dataOpac* dopa, const GLfloat *, GLfloat *buffer, int, bool=false )=0;
-    void     read_normals( int, int, const char * );
-    void     vertnorm( GLfloat *a ){ _ptnrml=a; }
-    int              ntris() { return _ptsPerObj-2; } //!< equivalent #triangles
-    
+    virtual int   buffer( int, GLfloat*, Colourscale*, DATA_TYPE*,
+                         dataOpac* dopa, const GLfloat *, GLfloat *buffer, int, char=0 )=0;
+    void          read_normals( int, int, const char * );
+    void          vertnorm( GLfloat *a ){ _ptnrml=a; }
+    int           ntris() { return _ptsPerObj-2; } //!< equivalent #triangles
   protected:
     GLfloat* _nrml;
     GLfloat* _ptnrml;
@@ -200,7 +204,7 @@ class PolyGon : public SurfaceElement
 
     virtual void compute_normals( int a, int b );
     virtual int  buffer( int, GLfloat*, Colourscale*, DATA_TYPE*,
-                         dataOpac* dopa, const GLfloat *, GLfloat*, int, bool=false );
+                         dataOpac* dopa, const GLfloat *, GLfloat*, int, char=0 );
     virtual bool read( const char * ){return true;}
 	static const int  _zero;
     const int* iso_polys(unsigned int index){return &_zero;}
@@ -214,7 +218,7 @@ class Triangle : public SurfaceElement
     Triangle(PPoint *p):SurfaceElement(p,3) {}
     virtual ~Triangle() {}
     virtual int      buffer( int, GLfloat*, Colourscale*, DATA_TYPE*,
-                         dataOpac*, const GLfloat *, GLfloat *, int, bool=false );
+                         dataOpac*, const GLfloat *, GLfloat *, int, char=0 );
     virtual bool     read( const char * );
     virtual void     compute_normals( int, int );
     virtual DrawingObj *isosurf( DATA_TYPE *d, DATA_TYPE val ){return NULL;}
@@ -233,7 +237,7 @@ class Quadrilateral : public SurfaceElement
     virtual ~Quadrilateral() {}
 
     virtual int      buffer( int, GLfloat*, Colourscale*, DATA_TYPE*,
-                         dataOpac*, const GLfloat *, GLfloat *, int, bool=false );
+                         dataOpac*, const GLfloat *, GLfloat *, int, char=0 );
     virtual bool     read( const char * );
     virtual void     compute_normals( int, int );
     bool     add( const char * );
diff --git a/src/Model.cc b/src/Model.cc
index 0331d06..84168f6 100644
--- a/src/Model.cc
+++ b/src/Model.cc
@@ -257,10 +257,7 @@ bool Model::read( const char* fnt, bool base1, bool no_elems )
   add_surface_from_surf( fn );
   LOG_TIMER("Add Surface from Surf");
 
-  //! \todo Figure out why we do this 
-  //for( int s=0; s<numSurf(); s++ )
-    //surface(s)->flip_norms();
-
+  associate_surf_with_vol();
   find_max_dim_and_bounds();
 
   return true;
@@ -841,6 +838,7 @@ int Model::add_surface_from_elem()
 	  continue;  //ignore volume elements
 	_surface[surfmap[reg]]->ele(surfs[reg])->define(idat);
 	_surface[surfmap[reg]]->ele(surfs[reg])->compute_normals(0,0);
+	_surface[surfmap[reg]]->src_ele(surfs[reg],i-1);
 	surfs[reg]++;
   }
 
@@ -1892,6 +1890,42 @@ Model::max_size( Object_t o )
 }
 
 
+/** determine the volume element to which a surface element is attached
+ */
+void
+Model::associate_surf_with_vol()
+{
+  map<int,set<int>> NtoV;  //!< map of volume elements associated with a node
+  for( int v=0; v<_vol.size(); v++ ) {
+    auto nodes = _vol[v]->obj();
+    for( int i = 0; i<_vol[v]->ptsPerObj(); i++ ) 
+      NtoV[nodes[i]].insert(v);
+  }
+
+  for( auto &surf : _surface ) {
+    auto elst = surf->ele();
+    for( int ei=0; ei<elst.size(); ei++ ) {
+      auto e=elst[ei];
+      if( surf->src_ele(ei) >= 0 ) continue;
+      int found;
+      for( auto n : NtoV[e->obj()[0]] ) {
+        found = n;
+        for( int i=1; i<e->ptsPerObj(); i++ ) { 
+            if( NtoV[e->obj()[i]].count(n) == 0 ) {
+              found = -1;
+              break;
+            }
+        }
+        if( found>=0 ) 
+          break;
+      }
+      surf->src_ele(ei,found);
+    }
+  }
+}
+
+
+
 #ifdef USE_HDF5
 /** read in one instant in time which is part of a element  time file
  *
diff --git a/src/Model.h b/src/Model.h
index 2f03974..24a887c 100644
--- a/src/Model.h
+++ b/src/Model.h
@@ -103,6 +103,7 @@ class Model
     bool             read_elem_file( const char* );
     void             increase_ele( int );
     bool             check_element( SurfaceElement *e );
+    void             associate_surf_with_vol();
 #ifdef USE_HDF5
     bool             read_elements(hid_t);
     bool             add_elements(hid_t hdf_file);
diff --git a/src/PolyGon.cc b/src/PolyGon.cc
index bd7c583..c85bc4c 100644
--- a/src/PolyGon.cc
+++ b/src/PolyGon.cc
@@ -12,14 +12,14 @@ const int PolyGon::_zero=0;
  *  \param nrml     if flat element normal, else vertex normals (NULL for none)
  *  \param buffer   buffer entry in which to place data
  *  \param surf     surface we are in 
- *  \param flat     flat shading
+ *  \param cmode    colouring mode
  *
  *  \return number of triangles put in draw buffer
  */
 int PolyGon::buffer( int p0, GLfloat *c, Colourscale* cs,
                           DATA_TYPE* data, dataOpac* dataopac,
                           const GLfloat* nrml, GLfloat *buffer, 
-                          int surf, bool flat )
+                          int surf, char cmode )
 {
   int i=p0*4;
   for ( int j=0; j<_ptsPerObj; j++ )
@@ -32,8 +32,10 @@ int PolyGon::buffer( int p0, GLfloat *c, Colourscale* cs,
 
       int k = j==2 ? 0 : n+1+j;
 
-      const GLfloat *col = colour( flat?p0:_node[i+k], c, data, cs, dataopac );  
-      buff_vtx( buff_off, _pt->pt(_node[i+k]), nrml+(flat?0:_node[i+k]*3), col, surf );
+      int cnode = cmode&(CM_ELEM|CM_FLAT) ? p0 : _node[i+k]; 
+      const GLfloat *col = colour( cnode, c, data, cs, dataopac );  
+      const GLfloat *n = nrml+((cmode&CM_FLAT)?0:_node[i+k]*3);
+      buff_vtx( buff_off, _pt->pt(_node[i+k]), n, col, surf );
     }
   return _ptsPerObj-2;
 }
diff --git a/src/Quadrilateral.cc b/src/Quadrilateral.cc
index 6aed4f4..6bd4e15 100644
--- a/src/Quadrilateral.cc
+++ b/src/Quadrilateral.cc
@@ -22,7 +22,7 @@ const int quad_iso_table[][11] = {
 
 /** draw one Quadrilaterals and do not do glBegin/glEnd
  *
- *  \param p0       index of first point to draw
+ *  \param p0       index of first quadrilateral to draw
  *  \param colour   colour to use if no data
  *  \param cs       colour scale
  *  \param data     data associated with nodes, element if flat (NULL for no data display)
@@ -30,12 +30,12 @@ const int quad_iso_table[][11] = {
  *  \param nrml     if flat, element normal, else vertex normals  (NULL for none)
  *  \param buffer   buffer entry in which to place data
  *  \param surf     surface we are in
- *  \param flat     flat shading
+ *  \param cmode    colour mode
  */
 int Quadrilateral::buffer( int p0, GLfloat *c, Colourscale* cs,
                           DATA_TYPE* data, dataOpac* dataopac,
                           const GLfloat* nrml,  GLfloat *buffer, 
-                                                  int surf, bool flat )
+                                              int surf, char cmode )
 {
   if ( p0>=_n ) return 0;
 
@@ -44,16 +44,20 @@ int Quadrilateral::buffer( int p0, GLfloat *c, Colourscale* cs,
        !_pt->vis(_node[i+2]) || !_pt->vis(_node[i+3])    )
     return 0;
 
-  GLfloat n[3] = { 0, 0, 0 };
+  GLfloat zn[3] = { 0, 0, 0 };
   for ( int j=0; j<3; j++ ) {
 
-    int node = _node[i+j];
-    const GLfloat *col = colour( flat?p0:node, c, data, cs, dataopac );  
-    buff_vtx( buffer+j*VBO_ELEM_SZ, _pt->pt(node), nrml?(flat?nrml:nrml+node*3):n, col, surf );
-
-    node = _node[i+(j+2)%4];
-    col =  colour( flat?p0:node, c, data, cs, dataopac );  
-    buff_vtx( buffer+(j+3)*VBO_ELEM_SZ, _pt->pt(node), nrml?(flat?nrml:nrml+node*3):n, col, surf );
+    int node  = _node[i+j];
+    int cnode = cmode&(CM_ELEM|CM_FLAT) ? p0 : node; 
+    const GLfloat *col = colour( cnode, c, data, cs, dataopac );  
+    const GLfloat *n = nrml?((cmode&CM_FLAT)?nrml:nrml+node*3):zn;
+    buff_vtx( buffer+j*VBO_ELEM_SZ, _pt->pt(node), n, col, surf );
+
+    node  = _node[i+(j+2)%4];
+    cnode = cmode&(CM_ELEM|CM_FLAT) ? p0 : node; 
+    col   = colour( cnode, c, data, cs, dataopac );  
+    n     = nrml?((cmode&CM_FLAT)?nrml:nrml+node*3):zn;
+    buff_vtx( buffer+(j+3)*VBO_ELEM_SZ, _pt->pt(node), n, col, surf );
   }
   return 2;
 }
diff --git a/src/Surfaces.cc b/src/Surfaces.cc
index 9296f32..52e8172 100644
--- a/src/Surfaces.cc
+++ b/src/Surfaces.cc
@@ -174,16 +174,16 @@ Surfaces::to_elem(GLfloat *ptdata)
  *  \param dat      if flat, element data, else nodal data (NULL for nodata display)
  *  \param dataopac data opacity
  *  \param buffer   buffer in which to place data
- *  \param flat     flat shading
+ *  \param cmode    colouring mode
  *
  *  \pre  the elements need to have been sorted by calling ::sort()
  *  \post the buffer reference points just past the entry added
  */
 void Surfaces::buffer( GLfloat *fill, Colourscale *cs, DATA_TYPE *dat,
-                     dataOpac* dataopac, GLfloat *&buffer, bool flat )
+                     dataOpac* dataopac, GLfloat *&buffer, char cmode )
 {
   for ( auto &a :_zlist )
-    buffer_elem( a.i, fill, cs, dat?(flat?dat+a.i:dat):NULL, dataopac, buffer, flat );
+    buffer_elem( a.i, fill, cs, dat, dataopac, buffer, cmode );
   _zlist.clear();
 }
 
@@ -196,13 +196,13 @@ void Surfaces::buffer( GLfloat *fill, Colourscale *cs, DATA_TYPE *dat,
  *  \param dat      nodal data (NULL for nodata display)
  *  \param dataopac data opacity
  *  \param buffer   buffer in which to place data
- *  \param flat     flat shading?
+ *  \param cmode    colouring mode
  *
  *  \post the buffer reference points just past the last entry added
  *  \post opaque elements are removed from \p _zlist
  */
 void Surfaces::buffer_opaque( float opaque, GLfloat *fill, Colourscale *cs, DATA_TYPE *dat,
-                     dataOpac* dataopac, GLfloat *&buffer, bool flat )
+                     dataOpac* dataopac, GLfloat *&buffer, char cmode )
 {
   set<int> solid;
 
@@ -218,7 +218,7 @@ void Surfaces::buffer_opaque( float opaque, GLfloat *fill, Colourscale *cs, DATA
       }
     if( eleopac ) {
       solid.insert(j);
-      buffer_elem( a.i, fill, cs, dat, dataopac, buffer, flat );
+      buffer_elem( a.i, fill, cs, dat, dataopac, buffer, cmode );
     }
   }
   for( auto i=solid.rbegin(); i != solid.rend(); i++ )
@@ -230,7 +230,7 @@ void Surfaces::buffer_opaque( float opaque, GLfloat *fill, Colourscale *cs, DATA
  * \param  ele      the element index
  *  \param fill     fill colour
  *  \param cs       colour scale
- *  \param dat      nodal data (NULL for nodata display)
+ *  \param dat      nodal or elemental data (NULL for nodata display)
  *  \param dataopac data opacity
  *  \param buf[out] buffer to fill
  *
@@ -238,10 +238,15 @@ void Surfaces::buffer_opaque( float opaque, GLfloat *fill, Colourscale *cs, DATA
  *  \post \p buf is advanced by the number of floats written
  */
 void Surfaces::buffer_elem( int ele, GLfloat *fill, Colourscale *cs, DATA_TYPE *dat,
-                            dataOpac* dataopac, GLfloat *&buf, bool flat )
+                            dataOpac* dataopac, GLfloat *&buf, char cmode )
 {
-  const GLfloat *nrml = flat?_ele[ele]->nrml():_vertnorm.get();
-  buf += _ele[ele]->buffer( 0, fill, cs, dat, dataopac, nrml, buf, _index, flat )*VBO_ELEM_SZ*3;
+  const GLfloat *nrml = (cmode&CM_FLAT)?_ele[ele]->nrml():_vertnorm.get();
+  GLfloat *dp = dat;
+  if( dat && cmode&CM_ELEM )
+    dp += _sele[ele];
+  else if( dat && cmode&CM_FLAT )
+    dp += ele;
+  buf += _ele[ele]->buffer( 0, fill, cs, dp, dataopac, nrml, buf, _index, cmode )*VBO_ELEM_SZ*3;
 }
 
 
diff --git a/src/Surfaces.h b/src/Surfaces.h
index 05d1c5f..9fbbba3 100644
--- a/src/Surfaces.h
+++ b/src/Surfaces.h
@@ -42,19 +42,19 @@ class Surfaces
     void get_vert_norms( GLfloat *&vn ){ vn = _vertnorm.get(); }
     void determine_vert_norms( PPoint & );
     SurfaceElement*& ele( int a ){ return _ele[a]; }
-    void   addele(int a,SurfaceElement*e){_ele[a]=e;}
+    void   addele(int a,SurfaceElement*e,int s=-1){_ele[a]=e;_sele[a]=s;}
     int    num() const {return _ele.size();}
-    void   num(int a){_ele.resize(a);}
+    void   num(int a){_ele.resize(a);_sele.resize(a, -1);}
     int    numtri(){ if(_numtri<0)count_tris();return _numtri; }  
     int    num_sharp_edges(){ return _sharp_edges.size(); }
     void   count_tris();
     vector<SurfaceElement*>& ele(){return _ele;}
     void   buffer( GLfloat *, Colourscale *, DATA_TYPE *,
-                     dataOpac*, GLfloat *&, bool=false );
+                     dataOpac*, GLfloat *&, char=0 );
     void   buffer_opaque( float, GLfloat *, Colourscale *, DATA_TYPE *,
-                     dataOpac*, GLfloat *&, bool = false );
+                     dataOpac*, GLfloat *&, char=0 );
     void   buffer_elem(int, GLfloat*,Colourscale*,DATA_TYPE*,
-                            dataOpac*, GLfloat *&, bool=false);
+                            dataOpac*, GLfloat *&, char=0);
     void   buffer_sharp_edges(GLfloat*&);
     void   zsort( void *, int, bool );
     void   label( string s ){ _label=s; }
@@ -71,16 +71,18 @@ class Surfaces
     GLfloat specular( void ){ return _specular[0]; }
     GLfloat shiny( void ) { return _shiny; }
     GLfloat backlight( void ) { return _backlight; }
-    void set_material( void );
+    void   set_material( void );
     vector<vtx_z>::iterator zl_begin() { return _zlist.begin(); }
     vector<vtx_z>::iterator zl_end() { return _zlist.end(); }
-    int  zl_sz(void){ return _zlist.size();}
-    void index( short i ){ _index = i; }
+    int    zl_sz(void){ return _zlist.size();}
+    void   index( short i ){ _index = i; }
     short  index(){ return _index; }
-    int  num_trans_tri() {int n=0; for(auto &l:_zlist)n+=_ele[l.i]->ntris(); return n;}
+    int     num_trans_tri() {int n=0; for(auto &l:_zlist)n+=_ele[l.i]->ntris(); return n;}
     vector<GLfloat> to_elem(GLfloat *ptdata);
     const unordered_set<unsigned int>& vertices(){return _vert;}
     pair<int,GLfloat*> sharp_interface( Surfaces * );
+    void    src_ele(int e, int s){ _sele[e] = s; }
+    long    src_ele(int a){ return _sele[a]; }
   protected:
     PPoint   *_p;
     GLfloat  _fillcolor[4]    = { 1., 0.5, 0.1, 1 };
@@ -91,6 +93,7 @@ class Surfaces
     unique_ptr<GLfloat[],free_delete> _vertnorm = NULL;   //!< vertex normals
     unordered_set<unsigned int>_vert;   //!< vertices for which normals are computed
     vector<SurfaceElement*> _ele;       //!< list of elelments to draw
+    vector<long> _sele;                 //!< source elements from which surface elements derived
     EdgeSet   _sharp_edges;             //!< sharp edges
     string   _label       = "";
     GLfloat  _diffuse[4]  = {0.6,0.6,0.6,1.};
diff --git a/src/TBmeshWin.cc b/src/TBmeshWin.cc
index 8a52dd8..44c24d0 100644
--- a/src/TBmeshWin.cc
+++ b/src/TBmeshWin.cc
@@ -27,6 +27,8 @@ const int MAX_SELECT_ATTEMPTS=3;
 #define REBUFFER_SURF(A) ( A & (VBO_Colour|VBO_Position|VBO_Visible) )
 #define REDRAW(A)  (_redraw_state&(A))
 
+#define CMODE(F,E) ((F)*CM_FLAT+(E)*CM_ELEM)
+
 unsigned int TBmeshWin::MAX_MESSAGES_READ = 100;
 
 
@@ -505,13 +507,14 @@ TBmeshWin::draw_surface(Surfaces* sf, GLfloat *&vbobuf, bool wf )
 
   glPolygonMode( GL_FRONT_AND_BACK, GL_FILL );
   sf->zsort( context(), stride, on_tr );
+  char cmode = CMODE( facetshading, dataBuffer->ele_based() );
   if( !on_tr ) 
     // all elements are opaque
-    sf->buffer( s_colour, cs, showData?(facetshading?sf->to_elem(data).data():data):NULL,
- 		NULL, vbobuf, facetshading );
+    sf->buffer( s_colour, cs, showData?((cmode==CM_FLAT)?sf->to_elem(data).data():data):NULL,
+ 		NULL, vbobuf, cmode );
   else if( translucency( obj, s_colour, true ) )
-    sf->buffer_opaque(  OPAQUE_LIMIT, s_colour, cs, facetshading?sf->to_elem(data).data():data,
-                                       dataopac->dop+obj, vbobuf, facetshading );
+    sf->buffer_opaque(  OPAQUE_LIMIT, s_colour, cs, (cmode==CM_FLAT)?sf->to_elem(data).data():data,
+                                       dataopac->dop+obj, vbobuf, cmode );
   if( !wf && on_tr ) {
     GLfloat *buff = _rd_shedge.buff_append( sf->num_sharp_edges()*2 );
     sf->buffer_sharp_edges( buff );
@@ -550,7 +553,7 @@ TBmeshWin::draw_sorted_elements( RenderTris &rd, vector<vtx_z> &elems, bool newT
           edata[e.s] = sf->to_elem(data);
       }
       sf->buffer_elem( e.i, scol, cs, showData?(facetshading?edata[e.s].data():data):NULL,
-              dataopac->dop+Surface, vbobuf, facetshading );
+              dataopac->dop+Surface, vbobuf, CMODE(facetshading,dataBuffer->ele_based()) );
     }
     rd.update_nodalbuff();
   }
@@ -1881,8 +1884,8 @@ TBmeshWin::draw_cut_planes( RRegion *reg )
       bool showData = cplane->datafied(i) && have_data!=NoData;
       clip_off( i );
       if( cut->fresh()                                ||                // first time VBO used after recalculated
-          _redraw_state&VBO_Clip                      ||                // clipping plane changed position
-          ( showData && _redraw_state&(VBO_Colour|VBO_Clip_Color)) ){             // data displayed on plane, and data changed
+              _redraw_state&VBO_Clip                      ||                // clipping plane changed position
+              ( showData && _redraw_state&(VBO_Colour|VBO_Clip_Color)) ){             // data displayed on plane, and data changed
         cut->init_buff(context());
         GLfloat *vbo = cut->allocate_buffer(); 
         for ( int e=0; e<cut->num(); e++ ) {
@@ -1892,18 +1895,25 @@ TBmeshWin::draw_cut_planes( RRegion *reg )
           for ( int v=1; v<cut->ele(e)->ptsPerObj(); v++ )
             memcpy( n+3*v, n, 3*sizeof(GLfloat ) );
 
-          // interpolate data
+          char cmode = CMODE( facetshading, dataBuffer->ele_based() );
+
           DATA_TYPE idata[cut->ele(e)->ptsPerObj()];
-          if ( showData && cplane->datafied(i) ) {
-            for ( int v=0; v<cut->ele(e)->ptsPerObj(); v++ ) {
-              if( _branch_cut )
-                idata[v] = cut->interpolate( e, data, v, _branch_range );
-              else
-                idata[v] = cut->interpolate( e, data, v );
+          if( cmode&CM_ELEM ) {
+            for ( int v=0; v<cut->ele(e)->ptsPerObj(); v++ ) 
+              idata[v] = data[cut->src_ele(e)];
+          }else {
+            // interpolate data
+            if ( showData && cplane->datafied(i) ) {
+              for ( int v=0; v<cut->ele(e)->ptsPerObj(); v++ ) {
+                if( _branch_cut )
+                  idata[v] = cut->interpolate( e, data, v, _branch_range );
+                else
+                  idata[v] = cut->interpolate( e, data, v );
+              }
             }
           }
           vbo += cut->ele(e)->buffer( 0, _cutplane_col, cs, showData?idata:NULL,
-                  dataopac->dop+Surface, n, vbo, 0 )*3*VBO_ELEM_SZ;
+                  dataopac->dop+Surface, n, vbo, 0, cmode )*3*VBO_ELEM_SZ;
         }
         cut->update_buff();
       } else if( _redraw_state&VBO_Clip_Color ) {           // update only opacity or undatified colour
@@ -1972,7 +1982,7 @@ TBmeshWin :: determine_cutplane( int cp )
         SurfaceElement *se = model->_vol[e]->cut( cpvis, cpf, interp );
         if ( se!= NULL ) {
 #pragma omp critical
-          _cutsurface[cp]->addEle( se, interp );
+          _cutsurface[cp]->addEle( se, interp, e );
         }
       }
     }
diff --git a/src/ThreadedData.h b/src/ThreadedData.h
index 4a35bf7..3e82ac4 100644
--- a/src/ThreadedData.h
+++ b/src/ThreadedData.h
@@ -149,7 +149,7 @@ class ThreadedData : public DataClass<T>
     virtual   T       max();	        //!< Maximum value of the entire series
     virtual   T       min(int);         //!< Minimum value at a slice of time
     virtual   T       min();            //!< Minimum value of the entire series
-    virtual   T*      slice(int);       //!< PPointer to slice of data at a time
+    virtual   T*      slice(int=-1);       //!< PPointer to slice of data at a time
     virtual   void    time_series( int, T* );    //!< PPointer to time series
     virtual   void    increment(int increment);  //!< Sets increment
   private:
diff --git a/src/Triangle.cc b/src/Triangle.cc
index 7d81662..e6b52b5 100644
--- a/src/Triangle.cc
+++ b/src/Triangle.cc
@@ -8,33 +8,35 @@ static const int tri_iso_table[][6] = {
 
 /** fill VBO with data for one Triangle
  *
- *  \param p0       index of first point to draw
+ *  \param p0       index of trianle to draw
  *  \param c        colour to use if no data
  *  \param cs       colour scale
- *  \param data     data associated with nodes, or elements if flat (NULL for no data display)
+ *  \param data     data associated with nodes or elements (NULL for no data display)
  *  \param dataopac data opacity (NULL for none)
  *  \param nrml     if \p flat, element normal, else vertex normals (NULL for none)
  *  \param buffer   buffer entry in which to place data
  *  \param surf     surface index
- *  \param flat     flat shading
+ *  \param cmode    colouring mode
  *  
  *  \return number of triangles buffered
  */
 int Triangle::buffer( int p0, GLfloat *c, Colourscale* cs,
                      DATA_TYPE* data, dataOpac* dataopac,
                      const GLfloat* nrml, GLfloat *buffer, 
-                                          int surf, bool flat )
+                                          int surf, char cmode )
 {
   if ( p0>=_n ) return 0;
   int i=3*p0;
   if ( !_pt->vis(_node[i]) || !_pt->vis(_node[i+1]) || !_pt->vis(_node[i+2]) )
     return 0;
 
-  GLfloat n[3] = { 0, 0, 0 };
+  GLfloat zn[3] = { 0, 0, 0 };
   for ( int j=0; j<3; j++ ) {
-    int nn = _node[i+j];
-    const GLfloat *col =  colour( flat?p0:nn, c, data, cs, dataopac );
-    buff_vtx( buffer+j*VBO_ELEM_SZ, _pt->pt(nn), nrml?(flat?nrml:nrml+nn*3):n, col, surf );
+    int            nn    = _node[i+j];
+    int            cnode = cmode&(CM_ELEM|CM_FLAT) ? p0 : nn; 
+    const GLfloat *col   =  colour( cnode, c, data, cs, dataopac );
+    const GLfloat *n     = nrml?((cmode&CM_FLAT)?nrml:nrml+nn*3):zn;
+    buff_vtx( buffer+j*VBO_ELEM_SZ, _pt->pt(nn), n, col, surf );
   }	    
   return 1;
 }
-- 
GitLab