[patch] Patch for the OpenFOAM Reader

OpenFOAMReader-3.2.1.patch (text/plain), 252.10 KB, created by Oliver Borm on 2007-12-18 12:13:56 UTC

--- ParaView3.orig/Servers/ServerManager/Resources/readers.xml	2007-10-20 15:05:53.381000000 +0900
+++ ParaView3/Servers/ServerManager/Resources/readers.xml	2007-10-20 13:31:17.732000000 +0900
@@ -3412,63 +3412,122 @@
    </SourceProxy>
 
     <!-- Beginning of FOAM Reader -->
-    <SourceProxy name="OpenFOAMReader" 
-                    class="vtkOpenFOAMReader">
-      <StringVectorProperty name="FileName"
-                            command="SetFileName"
-                            animateable="0"
-                            number_of_elements="1">
-        <FileListDomain name="files"/>
-        <Documentation>
-          Set the file name for the OpenFOAM reader.
-        </Documentation>
-      </StringVectorProperty>
-      <StringVectorProperty name="CellArrayInfo"
-                            information_only="1">
-        <ArraySelectionInformationHelper attribute_name="Cell"/>
-      </StringVectorProperty>
-      <StringVectorProperty name="CellArrayStatus"
-                            command="SetCellArrayStatus"
-                            number_of_elements="0" 
-                            repeat_command="1"
-                            number_of_elements_per_command="2"
-                            element_types="2 0" 
-                            information_property="CellArrayInfo"
-                            label="Cell Arrays">
-        <ArraySelectionDomain name="array_list">
-          <RequiredProperties>
-            <Property name="CellArrayInfo" function="ArrayList"/>
-          </RequiredProperties>
-        </ArraySelectionDomain>
-        <Documentation>
-          Select which cell-centered arrays to read.
-        </Documentation>
-      </StringVectorProperty>
-      <IntVectorProperty name="TimeStepRangeInfo"
-                         command="GetTimeStepRange"
-                         information_only="1">
-        <SimpleIntInformationHelper/>
-      </IntVectorProperty>
-      <IntVectorProperty name="TimeStep"
-                           command="SetTimeStep"
-                           number_of_elements="1"
-                           animateable="1"
-                           default_values="0">
-        <IntRangeDomain name="range">
-          <RequiredProperties>
-            <Property name="TimeStepRangeInfo" function="Range"/>
-          </RequiredProperties>
-        </IntRangeDomain>
-        <Documentation>
-          Set the current timestep.
-        </Documentation>
-      </IntVectorProperty>
-      <DoubleVectorProperty 
-          name="TimestepValues"
-          repeatable="1"
-          information_only="1">
-        <TimeStepsInformationHelper/>
-      </DoubleVectorProperty>
+    <SourceProxy
+
+        name="OpenFOAMReader" 
+        class="vtkOpenFOAMReader"
+	output="vtkDataSet">
+
+       <StringVectorProperty
+           name="FileName"
+           command="SetFileName"
+           number_of_elements="1">
+           <FileListDomain name="files"/>
+        </StringVectorProperty>
+
+	<IntVectorProperty 
+           name="CacheMesh" 
+           command="SetCacheMesh" 
+           number_of_elements="1"
+           default_values="1">
+           <BooleanDomain name="bool"/>
+                <Documentation>
+                   Cache the OpenFOAM mesh between GUI selection changes.
+                </Documentation>
+        </IntVectorProperty>
+
+        <IntVectorProperty name="ReadZones" 
+                           command="SetReadZones"
+                           number_of_elements="1" 
+                           default_values="0"
+                           label="Read Zones">
+          <BooleanDomain name="bool"/>
+          <Documentation>
+            Read point/face/cell-Zones?
+          </Documentation>
+        </IntVectorProperty>
+
+	    <IntVectorProperty name="AccumulatePatches" 
+			       command="SetKeepPatches"
+			       number_of_elements="1" 
+			       default_values="1"
+			       label="Accumulate Patches to Selection List">
+		    <BooleanDomain name="bool"/>
+		    <Documentation>
+			    Accumulate patches into selection list across time-steps even if they stop existing
+		    </Documentation>
+	    </IntVectorProperty>
+
+	<DoubleVectorProperty
+           name="TimestepValues"
+	   repeatable="1"
+	   information_only="1">
+           <TimeStepsInformationHelper/>
+        </DoubleVectorProperty>
+
+	<StringVectorProperty 
+            name="PatchArrayInfo"
+            information_only="1">
+            <ArraySelectionInformationHelper attribute_name="Patch"/>
+	</StringVectorProperty>
+	<StringVectorProperty
+            name="MeshRegions"
+            command="SetPatchArrayStatus"
+            number_of_elements="0" 
+            repeat_command="1"
+            number_of_elements_per_command="2"
+            element_types="2 0"
+            information_property="PatchArrayInfo">
+            <ArraySelectionDomain name="array_list">
+                <RequiredProperties>
+                    <Property name="PatchArrayInfo"
+                        function="ArrayList"/>
+                    </RequiredProperties>
+            </ArraySelectionDomain>
+        </StringVectorProperty>
+	
+        <StringVectorProperty 
+            name="CellArrayInfo"
+            information_only="1">
+            <ArraySelectionInformationHelper attribute_name="Cell"/>
+        </StringVectorProperty>
+        <StringVectorProperty
+            name="CellArrays"
+            command="SetCellArrayStatus"
+            number_of_elements="0" 
+            repeat_command="1"
+            number_of_elements_per_command="2"
+            element_types="2 0"
+    	    information_property="CellArrayInfo">
+            <ArraySelectionDomain name="array_list">
+                <RequiredProperties>
+                    <Property name="CellArrayInfo"
+                        function="ArrayList"/>
+                </RequiredProperties>
+            </ArraySelectionDomain>
+        </StringVectorProperty>
+
+        <StringVectorProperty 
+            name="PointArrayInfo"
+            information_only="1">
+            <ArraySelectionInformationHelper attribute_name="Point"/>
+        </StringVectorProperty>
+        <StringVectorProperty
+            name="PointArrays"
+            command="SetPointArrayStatus"
+            number_of_elements="0" 
+            repeat_command="1"
+            number_of_elements_per_command="2"
+            element_types="2 0"
+	        information_property="PointArrayInfo">
+            <ArraySelectionDomain name="array_list">
+                <RequiredProperties>
+                    <Property name="PointArrayInfo"
+                        function="ArrayList"/>
+                </RequiredProperties>
+            </ArraySelectionDomain>
+        </StringVectorProperty>
+
     </SourceProxy>
     <!-- End of foam Reader -->
 

--- ParaView3.orig/VTK/IO/vtkOpenFOAMReader.cxx	2007-11-07 21:51:56.000000000 +0100
+++ ParaView3/VTK/IO/vtkOpenFOAMReader.cxx	2007-10-23 12:37:05.000000000 +0200
@@ -16,11 +16,21 @@
 // Technology Laboratory who developed this class.
 // Please address all comments to Terry Jordan ([e-mail])
 //
+// Token-based FoamFile format lexer/parser, performance enhancements,
+// gzipped file and lagrangian field support by Takuya Oshima
+// ([e-mail])
+//
+// * GUI Based selection of mesh regions and fields available in the case
+// * Minor bug fixes / Strict memory allocation checks
+// * Minor performance enhancements
+// by Philippose Rajan ([e-mail])
+
 #include "vtkOpenFOAMReader.h"
 
+#include <vtksys/ios/sstream>
 #include <vtkstd/string>
 #include <vtkstd/vector>
-#include <vtksys/ios/sstream>
+#include <vtk_zlib.h>
 
 #include "vtkInformation.h"
 #include "vtkInformationVector.h"
@@ -28,10 +38,12 @@
 #include "vtkDataArraySelection.h"
 #include "vtkStreamingDemandDrivenPipeline.h"
 #include "vtkCellArray.h"
+#include "vtkCallbackCommand.h"
 #include "vtkDataArraySelection.h"
 #include "vtkIntArray.h"
-#include "vtkFloatArray.h"
 #include "vtkDoubleArray.h"
+#include "vtkStringArray.h"
+#include "vtkSortDataArray.h"
 #include "vtkPoints.h"
 #include "vtkCellData.h"
 #include "vtkHexahedron.h"
@@ -48,43 +60,2104 @@
 #include "vtkUnstructuredGridAlgorithm.h"
 #include "vtkObjectFactory.h"
 #include "vtkDirectory.h"
+#include "vtkPolyData.h"
+#include "vtkPointData.h"
+#include "vtkPVConfig.h" // for PARAVIEW_VERSION_MAJOR
 
 #include <ctype.h>
 #include <sys/stat.h>
 
+#if 0
 #ifdef VTK_USE_ANSI_STDLIB
 #define VTK_IOS_NOCREATE
 #else
 #define VTK_IOS_NOCREATE | ios::nocreate
 #endif
+#endif
+
+// The input buffer size in bytes. The purpose of the buffering is to
+// reduce the number of costly gzread() calls so the buffer size
+// doesn't have to be large.
+#define VTK_FOAMFILE_BUFSIZE (2048)
+
+// If the macro is set to 1 fast (but inaccurate) string to floating
+// point number conversion routine is used instead of the system
+// strtod(). This affects about 15% in performance for ascii cases.
+#define VTK_FOAMFILE_FAST_STRTOD 1
 
-vtkCxxRevisionMacro(vtkOpenFOAMReader, "$Revision: 1.10 $");
+vtkCxxRevisionMacro(vtkOpenFOAMReader, "$Revision: 1.7.2.2 $");
 vtkStandardNewMacro(vtkOpenFOAMReader);
 
-struct stdString
+struct stringVector
 {
-  vtkstd::string value;
+  vtkstd::vector< vtkStdString > value;
 };
 
-struct stringVector
+struct intVector
+{
+  vtkstd::vector< int > value;
+};
+
+struct doubleVector
+{
+  vtkstd::vector<double> value;
+};
+
+struct intVectorVector
+{
+  vtkstd::vector< vtkstd::vector< int > > value;
+};
+
+struct faceVectorVector
+{
+  vtkstd::vector< vtkstd::vector< face > > value;
+};
+
+struct unstructuredGridVector
+{
+  vtkstd::vector<vtkUnstructuredGrid*> value;
+};
+
+struct vtkOpenFOAMReader::vtkFoamError
+{
+private:
+  vtkStdString str_;
+public:
+  vtkFoamError(): str_() {}
+  vtkFoamError(vtkFoamError& e): str_(e.str()) {}
+  const vtkStdString str() const { return str_; }
+  // a super-easy way to make use of operator<<()'s defined in
+  // vtkstd::ostringstream class
+  template <class T> vtkFoamError& operator<<(const T& t)
+  { vtkstd::ostringstream os; os << t; str_ += os.str(); return *this; }
+};
+
+// token class
+// based on src/OpenFOAM/db/IOstreams/token/{token|tokenI}.H
+// a word token is treated as a string token.
+struct vtkOpenFOAMReader::vtkFoamToken
 {
-  vtkstd::vector< vtkstd::string > value;
+public:
+  enum tokenType { UNDEFINED, PUNCTUATION, STRING, LABEL, SCALAR, ERROR };
+private:
+  tokenType type_;
+
+  vtkFoamToken(const vtkFoamToken&);
+  void clear() { if(type_ == STRING) delete stringTokenPtr_; }
+
+public:
+  union
+  {
+    char punctuationToken_; vtkStdString* stringTokenPtr_; int labelToken_;
+    double scalarToken_;
+  };
+
+  vtkFoamToken(): type_(UNDEFINED) {}
+  ~vtkFoamToken() { clear(); }
+  void setBad() { clear(); type_ = ERROR; }
+  tokenType type() const { return type_; }
+  bool isNumber() const { return type_ == LABEL || type_ == SCALAR; }
+
+  char punctuationToken() const { return punctuationToken_; }
+  int labelToken() const { return labelToken_; }
+  double number() const
+  { if(type_ == LABEL) return labelToken_; else return scalarToken_; }
+  const vtkStdString& stringToken() const
+  { return *stringTokenPtr_; }
+
+  void operator=(const char c)
+  { clear(); type_ = PUNCTUATION; punctuationToken_ = c; }
+  void operator=(const char *s)
+  { clear(); type_ = STRING; stringTokenPtr_ = new vtkStdString(s); }
+  void operator=(const vtkStdString& s)
+  { clear(); type_ = STRING; stringTokenPtr_ = new vtkStdString(s); }
+  void operator=(const int l)
+  { clear(); type_ = LABEL; labelToken_ = l; }
+  void operator=(const double d)
+  { clear(); type_ = SCALAR; scalarToken_ = d; }
+  bool operator==(const char c) const
+  { return type_ == PUNCTUATION && punctuationToken_ == c; }
+  bool operator==(const int i) const
+  { return type_ == LABEL && labelToken_ == i; }
+  bool operator==(const vtkStdString& str) const
+  { return type_ == STRING && *stringTokenPtr_ == str; }
+  bool operator!=(const vtkStdString& str) const
+  { return type_ != STRING || *stringTokenPtr_ != str; }
+  bool operator!=(const char c) const { return !operator==(c); }
+
+  friend vtkstd::ostringstream& operator<<(vtkstd::ostringstream& str,
+    const vtkFoamToken& t)
+  {
+    if(t.type() == PUNCTUATION)
+      {
+      str << t.punctuationToken_;
+      }
+    else if(t.type() == STRING)
+      {
+      str << *t.stringTokenPtr_;
+      }
+    else if(t.type() == LABEL)
+      {
+      str << t.labelToken_;
+      }
+    else if(t.type() == SCALAR)
+      {
+      str << t.scalarToken_;
+      }
+    else if(t.type() == ERROR)
+      {
+      str << "badToken (an unexpected EOF?)";
+      }
+    return str;
+  }
 };
 
-struct intVector
+// read and tokenize the input
+// line number counting is based on
+// src/OpenFOAM/db/IOstreams/Sstreams/ISstreamI.H
+struct vtkOpenFOAMReader::vtkFoamFile
+{
+private:
+  gzFile file_;
+  bool putBack_;
+  int putBackC_;
+  int lineNumber_;
+
+  // buffer pointers. using raw pointers for performance reason.
+  unsigned char *buf_;
+  unsigned char *bufPtr_;
+  unsigned char *bufEndPtr_;
+
+  // declare and define as private
+  void putBack(const int c)
+  {
+    if(putBack_)
+      {
+      throw vtkFoamError() << "Attempted duplicated putBack()";
+      }
+    else
+      {
+      if(c == '\n')
+        {
+        lineNumber_--;
+        }
+      *--bufPtr_ = c;
+      putBack_ = true;
+      }
+  }
+
+  // get a character
+  int getc()
+  {
+    putBack_ = false;
+    if(bufPtr_ == bufEndPtr_)
+      {
+      bufPtr_ = buf_ + 1; // reserve the first byte for getback char
+      const int readSize = gzread(file_, bufPtr_, VTK_FOAMFILE_BUFSIZE - 1);
+      if(readSize <= 0) // 0 byte read or EOF
+        {
+        // set the current location to the end of the buffer so that
+        // getc() returns EOF again when called next time
+        bufPtr_ = bufEndPtr_;
+        return EOF;
+        }
+      bufEndPtr_ = bufPtr_ + readSize;
+      }
+    const int c = *bufPtr_++;
+    if(c == '\n')
+      {
+      lineNumber_++;
+      }
+    return c;
+  }
+
+  // get next semantically valid character
+  // based on src/OpenFOAM/db/IOstreams/Sstream/ISnextValid.C
+  int nextValid()
+  {
+    int c;
+    for(;;)
+      {
+      if((c = getc()) != EOF && isspace(c))
+        {
+        while((c = getc()) != EOF && isspace(c));
+        }
+      if(c == '/')
+        {
+        if((c = getc()) == EOF)
+          {
+          return '/';
+          }
+        if(c == '/') // one-line comment
+          {
+          while(EOF != (c = getc()) && c != '\n' && c != '\r');
+          }
+        else if(c == '*') // multi-line comment
+          {
+          for(;;)
+            {
+            if((c = getc()) == '*')
+              {
+              if((c = getc()) == '/')
+                {
+                break;
+                }
+              else if(c == EOF)
+                {
+                return EOF;
+                }
+              else
+                {
+                putBack(c);
+                }
+              }
+            else if(c == EOF)
+              {
+              return EOF;
+              }
+            }
+          }
+        else
+          {
+          putBack(c);
+          return '/';
+          }
+        }
+      else // a valid character or an EOF
+        {
+        return c;
+        }
+      }
+  }
+
+  // valid as a character that constitutes a word or not
+  // based on src/OpenFOAM/primitives/strings/word/wordI.H
+  bool valid(const int c) const
+  {
+    return (!isspace(c) && c != '"' && c != '/' && c != ';' && c != '{'
+      && c != '}');
+  }
+
+public:
+  vtkFoamFile(): file_(NULL), putBack_(false), putBackC_(0), lineNumber_(0),
+    buf_(NULL) {}
+  ~vtkFoamFile() { close(); }
+
+  void open(const vtkStdString& path)
+  {
+    lineNumber_ = 0;
+
+    buf_ = new unsigned char[VTK_FOAMFILE_BUFSIZE];
+    bufPtr_ = buf_;
+    bufEndPtr_ = buf_;
+
+    if(file_ != NULL)
+      {
+      throw vtkFoamError() << "File already opened";
+      }
+    if((file_ = gzopen(path.c_str(), "rb")) == NULL)
+      {
+      throw vtkFoamError() << "Can't open";
+      }
+
+    lineNumber_ = 1;
+  }
+
+  void close()
+  {
+    if(buf_ != NULL)
+      {
+      delete [] buf_;
+      buf_ = NULL;
+      }
+    putBack_ = false;
+    putBackC_ = 0;
+    if(file_ != NULL)
+      {
+      gzclose(file_);
+      file_ = NULL;
+      }
+    // don't reset the line number so that the last line number is
+    // retained after close
+    // lineNumber_ = 0;
+  }
+
+  int lineNumber() const { return lineNumber_; }
+
+  // gzread with buffering handling
+  int read(char *buf, const int len)
+  {
+    int readlen;
+    const int buflen = bufEndPtr_ - bufPtr_;
+    if(len > buflen)
+      {
+      memcpy(buf, bufPtr_, buflen);
+      readlen = gzread(file_, buf + buflen, len - buflen);
+      if(readlen >= 0)
+        {
+        readlen += buflen;
+        }
+      else
+        {
+        if(buflen == 0) // return EOF
+          {
+          readlen = -1;
+          }
+        else
+          {
+          readlen = buflen;
+          }
+        }
+      bufPtr_ = bufEndPtr_;
+      }
+    else
+      {
+      memcpy(buf, bufPtr_, len);
+      bufPtr_ += len;
+      readlen = len;
+      }
+    for(int i = 0; i < readlen; i++)
+      {
+      if(buf[i] == '\n')
+        {
+        lineNumber_++;
+        }
+      }
+    return readlen;
+  }
+
+  // the tokenizer
+  // based on src/OpenFOAM/db/IOstreams/Sstreams/{ISreadToken|ISread}.C
+  // returns true if success, false if encountered EOF
+  bool read(vtkFoamToken& token)
+  {
+    int c = nextValid();
+    if(c == EOF)
+      {
+      token.setBad();
+      return false;
+      }
+    switch(c)
+      {
+      case ';': case '(': case ')': case '{': case '}': case '[': case ']':
+      case ':': case ',': case '=': case '+': case '*': case '/':
+        // punctuation token
+        token = (char)c;
+        return true;
+      case '-': case '.': case '0': case '1': case '2': case '3': case '4':
+      case '5': case '6': case '7': case '8': case '9':
+        // number token
+        {
+        const int MAX_NUM = 100;
+        static char numberBuffer[MAX_NUM];
+        bool isScalar = false;
+        if(c == '.')
+          {
+          isScalar = true;
+          }
+        int i = 0, isDigit;
+        numberBuffer[i++] = c;
+        while(EOF != (c = getc()) && i < MAX_NUM && ((isDigit = isdigit(c))
+          || c == '.' || c == 'e' || c == 'E' || c == '+' || c == '-'))
+          {
+          numberBuffer[i++] = c;
+          if(!isDigit)
+            {
+            isScalar = true;
+            }
+          }
+        if(i == MAX_NUM)
+          {
+          if(c != EOF)
+            {
+            putBack(c);
+            }
+          numberBuffer[MAX_NUM - 1] = '\0';
+          token.setBad();
+          throw vtkFoamError() << "Reached the maximum allowed number length";
+          }
+        numberBuffer[i] = '\0';
+        if(c != EOF)
+          {
+          putBack(c);
+          if(i == 1 && numberBuffer[0] == '-')
+            {
+            token = '-';
+            }
+          else if(isScalar)
+            {
+            token = strtod(numberBuffer, NULL);
+            }
+          else
+            {
+            token = static_cast<int>(strtol(numberBuffer, NULL, 10));
+            }
+          return true;
+          }
+        else
+          {
+          throw vtkFoamError() << "Encountered EOF while tokenizing a number";
+          }
+        }
+      case '"':
+        // string token
+        {
+        const int MAX_STR = 1024;
+        static char stringBuffer[MAX_STR];
+        int i = 0;
+        bool escape = false;
+        while((c = getc()) != EOF)
+          {
+          if(c == '"' && !escape)
+            {
+            stringBuffer[i] = '\0';
+            token = stringBuffer;
+            return true;
+            }
+          if((c == '\n' || c == '\r') && !escape)
+            {
+            stringBuffer[i] = '\0';
+            token = stringBuffer;
+            throw vtkFoamError() << "Found a newline while reading string \""
+                << stringBuffer << "\"";
+            }
+          if(c == '\\')
+            {
+            escape = true;
+            }
+          else
+            {
+            escape = false;
+            stringBuffer[i] = c;
+            if(++i == MAX_STR)
+              {
+              stringBuffer[MAX_STR - 1] = '\0';
+              token = stringBuffer;
+              throw vtkFoamError()
+                << "Reached the maximum allowed string length";
+              }
+            }
+          }
+        stringBuffer[i] = '\0';
+        token = stringBuffer;
+        throw vtkFoamError() << "Encountered EOF while reading string";
+        }
+      default:
+        // parses as a word token, but gives the STRING type for simplicity
+        {
+        putBack(c);
+        const int MAX_WORD = 1024;
+        static char wordBuffer[MAX_WORD];
+        int i = 0;
+        int bc = 0;
+        while((EOF != (c = getc())) && valid(c))
+          {
+          if(i >= MAX_WORD)
+            {
+            wordBuffer[MAX_WORD - 1] = '\0';
+            token = wordBuffer;
+            throw vtkFoamError() << "Reached the maximum allowed word length.";
+            }
+          if(c == '(')
+            {
+            bc++;
+            }
+          else if(c == ')')
+            {
+            bc--;
+            if(bc == -1)
+              {
+              break;
+              }
+            }
+          wordBuffer[i++] = c;
+          }
+        if(i == 0)
+          {
+          token.setBad();
+          throw vtkFoamError() << "Invalid first character as a token";
+          }
+        wordBuffer[i] = '\0';
+        token = wordBuffer;
+        if(c != EOF)
+          {
+          putBack(c);
+          }
+        return true;
+        }
+      }
+  }
+
+  void readExpecting(const char c)
+  {
+    vtkFoamToken t;
+    if(!read(t))
+      {
+      throw vtkFoamError() << "Unexpected EOF";
+      }
+    if(t != c)
+      {
+      throw vtkFoamError() << "Expected punctuation token " << c << ", found "
+        << t;
+      }
+  }
+
+  void readExpecting(const char* str)
+  {
+    vtkFoamToken t;
+    if(!read(t))
+      {
+      throw vtkFoamError() << "Unexpected EOF";
+      }
+    if(t != str)
+      {
+      throw vtkFoamError() << "Expected string \"" << str << "\", found " << t;
+      }
+  }
+
+  // specialized for reading an integer value.
+  // not using the standard strtol() for speed reason.
+  int readLabel()
+  {
+    int c = nextValid();
+
+    bool negative = false;
+    if(c == '-')
+      {
+      negative = true;
+      c = getc();
+      }
+
+    if(!isdigit(c))
+      {
+      throw vtkFoamError()
+        << "Expected an integer, found a non-digit character" << (char)c;
+      }
+
+    int num = c - '0';
+    while(EOF != (c = getc()) && isdigit(c))
+      {
+      num = 10 * num + c - '0';
+      }
+
+    if(c == EOF)
+      {
+      throw vtkFoamError() << "Encountered EOF while reading integer";
+      }
+    putBack(c);
+
+    return negative ? -num : num;
+  }
+
+#if VTK_FOAMFILE_FAST_STRTOD
+  // extreamely simplified high-performing string to floating point
+  // conversion code based on
+  // ParaView3/VTK/Utilities/vtksqlite/vtk_sqlite3.c
+  double readNumber()
+  {
+    int c = nextValid();
+
+    // determine sign
+    bool negative = false;
+    if(c == '-')
+      {
+      negative = true;
+      c = getc();
+      }
+
+    if(c == EOF || (!isdigit(c) && c != '.' && c != 'e' && c != 'E'))
+      {
+      throw vtkFoamError()
+        << "Expected a number, found a non-digit character" << (char)c;
+      }
+
+    // read integer part
+    long double num = 0.0;
+    while(c != EOF && isdigit(c))
+      {
+      num = num * 10.0 + (c - '0');
+      c = getc();
+      }
+
+    // read decimal part
+    if(c == '.')
+      {
+      long double divisor = 1.0;
+
+      c = getc();
+      while(c != EOF && isdigit(c))
+        {
+        num = num * 10.0 + (c - '0');
+        divisor *= 10.0;
+        c = getc();
+        }
+      num /= divisor;
+      }
+
+    // read exponent part
+    if(c == 'e' || c == 'E')
+      {
+      int esign = 1;
+      int eval = 0;
+      long double scale = 1.0;
+
+      c = getc();
+      if(c == '-')
+        {
+        esign = -1;
+        c = getc();
+        }
+      else if(c == '+')
+        {
+        c = getc();
+        }
+
+      while(c != EOF && isdigit(c))
+        {
+        eval = eval * 10 + c - '0';
+        c = getc();
+        }
+
+      // fast exponent multiplication!
+      while(eval >= 64)
+        {
+        scale *= 1.0e+64;
+        eval -= 64;
+        }
+      while(eval >= 16)
+        {
+        scale *= 1.0e+16;
+        eval -= 16;
+        }
+      while(eval >= 4)
+        {
+        scale *= 1.0e+4;
+        eval -= 4;
+        }
+      while(eval >= 1)
+        {
+        scale *= 1.0e+1;
+        eval -= 1;
+        }
+
+      if(esign < 0)
+        {
+        num /= scale;
+        }
+      else
+        {
+        num *= scale;
+        }
+      }
+
+    if(c == EOF)
+      {
+      throw vtkFoamError() << "Unexpected EOF";
+      }
+    putBack(c);
+
+    return static_cast<double>(negative ? -num : num);
+  }
+#else
+  // specialized for reading a scalar value
+  double readNumber()
+  {
+    int c = nextValid();
+    if(c == EOF || (!isdigit(c) && c != '-' && c != '.'))
+      {
+      if(c == EOF)
+        {
+        throw vtkFoamError() << "Unexpected EOF";
+        }
+      else
+        {
+        throw vtkFoamError() << "Expected a number, found " << (char)c;
+        }
+      }
+
+    const int MAX_NUM = 100;
+    static char numberBuffer[MAX_NUM];
+    int i = 0;
+    numberBuffer[i++] = c;
+    if(EOF != (c = getc()) && i < MAX_NUM && (isdigit(c)
+      || c == '.' || c == 'e' || c == 'E')) // signs do not come as a 2nd char
+      {
+      numberBuffer[i++] = c;
+      while(EOF != (c = getc()) && i < MAX_NUM && (isdigit(c)
+        || c == '.' || c == 'e' || c == 'E' || c == '+' || c == '-'))
+        {
+        numberBuffer[i++] = c;
+        }
+      }
+    if(i == MAX_NUM)
+      {
+      throw vtkFoamError() << "Reached the maximum allowed number length";
+      }
+    numberBuffer[i] = '\0';
+    if(c == EOF)
+      {
+      throw vtkFoamError() << "Unexpected EOF";
+      }
+    else
+      {
+      putBack(c);
+      if(i == 1 && numberBuffer[0] == '-')
+        {
+        throw vtkFoamError() << "Expected a number, found a minus sign";
+        }
+      }
+    return strtod(numberBuffer, NULL);
+  }
+#endif
+};
+
+// IOobject class
+// holds file handle, file format, name of the object the file holds and
+// type of the object.
+struct vtkOpenFOAMReader::vtkFoamIOobject: public vtkFoamFile
+{
+public:
+  enum fileFormat { UNDEFINED, ASCII, BINARY };
+
+private:
+  vtkStdString fileName_;
+  fileFormat format_;
+  vtkStdString objectName_;
+  vtkStdString headerClassName_;
+  vtkFoamError* e_;
+
+  void readHeader(); // defined later
+public:
+  vtkFoamIOobject(): vtkFoamFile(), format_(UNDEFINED), e_(NULL) {}
+  ~vtkFoamIOobject() { close(); }
+
+  bool open(const vtkStdString& file)
+  {
+    fileName_ = file;
+    try
+      {
+      vtkFoamFile::open(file);
+      }
+    catch(vtkFoamError& e)
+      {
+      e_ = new vtkFoamError(e);
+      return false;
+      }
+
+    try
+      {
+      readHeader();
+      }
+    catch(vtkFoamError& e)
+      {
+      vtkFoamFile::close();
+      e_ = new vtkFoamError(e);
+      return false;
+      }
+    e_ = NULL;
+    return true;
+  }
+
+  void close()
+  {
+    vtkFoamFile::close();
+    format_ = UNDEFINED;
+    objectName_.erase();
+    headerClassName_.erase();
+    if(e_ != NULL)
+      {
+      delete e_;
+      }
+  }
+  const vtkStdString& fileName() { return fileName_; }
+  const fileFormat format() const { return format_; }
+  const vtkStdString& className() const { return headerClassName_; }
+  const vtkStdString& objectName() const { return objectName_; }
+  const vtkFoamError& error() const { return *e_; }
+  void setError(vtkFoamError& e)
+  { if(e_ != NULL) delete e_; e_ = new vtkFoamError(e); }
+};
+
+// a class that represents a value of a dictionary entry that corresponds to
+// its keyword. note that an entry can have more than one value.
+struct vtkOpenFOAMReader::vtkFoamEntryValue
+{
+public:
+  enum entryValueType
+    {
+    UNDEFINED, PUNCTUATION, STRING, LABEL, SCALAR, STRINGLIST, LABELLIST,
+    LABELLISTLIST, SCALARLIST, VECTORLIST, ENTRYVALUELIST, EMPTYLIST,
+    DICTIONARY, ERROR
+    };
+
+private:
+  entryValueType type_;
+  bool isUniform_;
+  bool managed_;
+
+  // don't pack them into a union
+  vtkFoamToken t_;
+  vtkstd::vector<vtkFoamEntryValue*> entryValuePtrs_;
+
+  union
+  {
+    intVector *labelListPtr_; vtkDoubleArray *scalarListPtr_, *vectorListPtr_;
+    stringVector *stringListPtr_; intVectorVector *labelListListPtr_;
+    vtkFoamDict *dictPtr_;
+  };
+
+  void clear();
+  void readList(vtkFoamIOobject& io);
+
+public:
+  vtkFoamEntryValue(): type_(UNDEFINED), isUniform_(false), managed_(true),
+    labelListPtr_(NULL) {}
+  ~vtkFoamEntryValue() { clear(); }
+
+  entryValueType type() const { return type_; }
+  void setEmptyList() { clear(); isUniform_ = false; type_ = EMPTYLIST; }
+  bool isUniform() const { return isUniform_; }
+  void read(vtkFoamIOobject& io);
+  void readDictionary(vtkFoamIOobject& io,
+    const vtkStdString& firstKeyword = "");
+  const vtkFoamToken& token() const { return t_; }
+  vtkFoamToken& token() { return t_; }
+  const vtkstd::vector<int>& labelList() const
+  { return labelListPtr_->value; }
+  const vtkstd::vector<vtkstd::vector<int> >& labelListList() const
+  { return labelListListPtr_->value; }
+  vtkDoubleArray& scalarList()
+  { return *scalarListPtr_; }
+  vtkDoubleArray& vectorList()
+  { return *vectorListPtr_; }
+  vtkFoamDict& dictionary()
+  { return *dictPtr_; }
+  bool operator==(const char c) const
+  { return type_ == PUNCTUATION && t_ == c; }
+  bool operator!=(const char c) const
+  { return type_ != PUNCTUATION || t_ != c; }
+  bool operator==(const int i) const
+  { return type_ == LABEL && t_ == i; }
+
+  void *ptr()
+  {
+    managed_ = false; // the returned pointer will not be deleted by the d'tor
+    return (void *)labelListPtr_; // all list pointers are in a single union
+  }
+
+  // the following two are for an exceptional expression of
+  // `LABEL{LABELorSCALAR}' without type prefix (e. g. `2{-0}' in
+  // mixedRhoE B.C. in rhopSonicFoam/shockTube)
+  void makeLabelList(const int labelValue, const int size)
+  {
+    labelListPtr_ = new intVector;
+    type_ = LABELLIST;
+    vtkstd::vector<int>& ll = labelListPtr_->value;
+    ll.resize(size);
+    for(int i = 0; i < size; i++)
+      {
+      ll[i] = labelValue;
+      }
+  }
+  void makeScalarList(const double scalarValue, const int size)
+  {
+    scalarListPtr_ = vtkDoubleArray::New();
+    type_ = SCALARLIST;
+    vtkDoubleArray& sl = *scalarListPtr_;
+    sl.SetNumberOfValues(size);
+    for(vtkIdType i = 0; i < size; i++)
+      {
+      sl.SetValue(i, scalarValue);
+      }
+  }
+
+  // read a labelList
+  void readLabelList(vtkFoamIOobject& io)
+  {
+    vtkFoamToken currToken;
+    if(!io.read(currToken))
+      {
+      throw vtkFoamError() << "Unexpected EOF";
+      }
+    labelListPtr_ = new intVector;
+    type_ = LABELLIST;
+    if(currToken.type() == vtkFoamToken::LABEL)
+      {
+      const int size = currToken.labelToken();
+      labelListPtr_->value.resize(size);
+      if(!io.read(currToken))
+        {
+        throw vtkFoamError() << "Unexpected EOF";
+        }
+      vtkstd::vector<int>& ll = labelListPtr_->value;
+      // some objects have lists with only one element enclosed by {}
+      // e. g. simpleFoam/pitzDaily3Blocks/constant/polyMesh/faceZones
+      if(currToken == '{')
+        {
+        int labelValue = io.readLabel();
+        for(int i = 0; i < size; i++)
+          {
+          ll[i] = labelValue;
+          }
+        io.readExpecting('}');
+        return;
+        }
+      if(currToken != '(')
+        {
+        throw vtkFoamError() << "Expected (, found " << currToken;
+        }
+      if(io.format() == vtkFoamIOobject::ASCII)
+        {
+        for(int i = 0; i < size; i++)
+          {
+          ll[i] = io.readLabel();
+          }
+        }
+      else
+        {
+        if(size > 0) // avoid invalid access to ll.at(0)
+          {
+          io.read(reinterpret_cast<char*>(&ll.at(0)), size * sizeof(int));
+          }
+        }
+      io.readExpecting(')');
+      }
+    else if(currToken == '(')
+      {
+      while(io.read(currToken) && currToken != ')')
+        {
+        if(currToken.type() != vtkFoamToken::LABEL)
+          {
+          throw vtkFoamError() << "Expected an integer or a (, found "
+            << currToken;
+          }
+        labelListPtr_->value.push_back(currToken.labelToken());
+        }
+      }
+    else
+      {
+      throw vtkFoamError() << "Expected integer or (, found " << currToken;
+      }
+  }
+
+  // reads a list of labelLists. requires size prefix of the listList
+  // to be present. size of the list must be present in the stream if
+  // the format is binary.
+  void readLabelListList(vtkFoamIOobject& io)
+  {
+    vtkFoamToken currToken;
+    if(!io.read(currToken))
+      {
+      throw vtkFoamError() << "Unexpected EOF";
+      }
+    if(currToken.type() == vtkFoamToken::LABEL)
+      {
+      labelListListPtr_ = new intVectorVector;
+      type_ = LABELLISTLIST;
+
+      const int sizeI = currToken.labelToken();
+      labelListListPtr_->value.resize(sizeI);
+      io.readExpecting('(');
+      for(int i = 0; i < sizeI; i++)
+        {
+        if(!io.read(currToken))
+          {
+          throw vtkFoamError() << "Unexpected EOF";
+          }
+        if(currToken.type() == vtkFoamToken::LABEL)
+          {
+          const int sizeJ = currToken.labelToken();
+          labelListListPtr_->value[i].resize(sizeJ);
+          io.readExpecting('(');
+          vtkstd::vector<int>& labelListI = labelListListPtr_->value[i];
+          if(io.format() == vtkFoamIOobject::ASCII)
+            {
+            for(int j = 0; j < sizeJ; j++)
+              {
+              labelListI[j] = io.readLabel();
+              }
+            }
+          else
+            {
+            if(sizeJ > 0) // avoid invalid reference to labelListI.at(0)
+              {
+              io.read(reinterpret_cast<char*>(&labelListI.at(0)),
+                sizeJ * sizeof(int));
+              }
+            }
+          io.readExpecting(')');
+          }
+        else if(currToken == '(')
+          {
+          while(io.read(currToken) && currToken != ')')
+            {
+            if(currToken.type() != vtkFoamToken::LABEL)
+              {
+              throw vtkFoamError() << "Expected an integer, found "
+                << currToken;
+              }
+            labelListListPtr_->value[i].push_back(currToken.labelToken());
+            }
+          }
+        else
+          {
+          throw vtkFoamError() << "Expected integer or (, found " << currToken;
+          }
+        }
+      io.readExpecting(')');
+      }
+    else
+      {
+      throw vtkFoamError() << "Expected integer, found " << currToken;
+      }
+  }
+
+  // reads a scalarList. requires size prefix of the list to be
+  // present in the stream if the format is binary.
+  void readScalarList(vtkFoamIOobject& io)
+  {
+    vtkFoamToken currToken;
+    if(!io.read(currToken))
+      {
+      throw vtkFoamError() << "Unexpected EOF";
+      }
+    scalarListPtr_ = vtkDoubleArray::New();
+    type_ = SCALARLIST;
+    if(currToken.type() == vtkFoamToken::LABEL)
+      {
+      const vtkIdType size = currToken.labelToken();
+      scalarListPtr_->SetNumberOfValues(size);
+      if(!io.read(currToken))
+        {
+        throw vtkFoamError() << "Unexpected EOF";
+        }
+      vtkDoubleArray& sl = *scalarListPtr_;
+      // some objects have lists with only one element enclosed by {}
+      // to represent a field with an uniform value
+      if(currToken == '{')
+        {
+        double scalarValue = io.readNumber();
+        for(vtkIdType i = 0; i < size; i++)
+          {
+          sl.SetValue(i, scalarValue);
+          }
+        io.readExpecting('}');
+        return;
+        }
+      if(currToken != '(')
+        {
+        throw vtkFoamError() << "Expected (, found " << currToken;
+        }
+      if(io.format() == vtkFoamIOobject::ASCII)
+        {
+        for(vtkIdType i = 0; i < size; i++)
+          {
+          sl.SetValue(i, io.readNumber());
+          }
+        }
+      else
+        {
+        if(size > 0) // avoid invalid access to GetPointer(0)
+          {
+          io.read(reinterpret_cast<char*>(sl.GetPointer(0)),
+            size * sizeof(double));
+          }
+        }
+      io.readExpecting(')');
+      }
+    else if(currToken == '(')
+      {
+      while(io.read(currToken) && currToken != ')')
+        {
+        if(!io.read(currToken) || !currToken.isNumber())
+          {
+          throw vtkFoamError() << "Expected an integer or a (, found "
+            << currToken;
+          }
+        scalarListPtr_->InsertNextValue(currToken.number());
+        }
+      scalarListPtr_->Squeeze();
+      }
+    else
+      {
+      throw vtkFoamError() << "Expected integer or (, found " << currToken;
+      }
+  }
+
+  // reads a vectorList. requires size prefix of the list to be
+  // present in the stream if the format is binary.
+  // if isPositions is true read Cloud type of data as particle positions
+  // based on src/lagrangian/basic/particle/particleIO.C
+  void readVectorList(vtkFoamIOobject& io, const bool isPositions)
+  {
+    vtkFoamToken currToken;
+    if(!io.read(currToken))
+      {
+      throw vtkFoamError() << "Unexpected EOF";
+      }
+    vectorListPtr_ = vtkDoubleArray::New();
+    vectorListPtr_->SetNumberOfComponents(3);
+    type_ = VECTORLIST;
+    if(currToken.type() == vtkFoamToken::LABEL)
+      {
+      const int size = currToken.labelToken();
+      vectorListPtr_->SetNumberOfTuples(size);
+      if(!io.read(currToken))
+        {
+        throw vtkFoamError() << "Unexpected EOF";
+        }
+      // some objects may have a vector list with one element enclosed by {}
+      if(currToken == '{')
+        {
+        io.readExpecting('(');
+        double vectorValue[3];
+        vectorValue[0] = io.readNumber();
+        vectorValue[1] = io.readNumber();
+        vectorValue[2] = io.readNumber();
+        for(int i = 0; i < size; i++)
+          {
+          vectorListPtr_->SetTuple(i, vectorValue);
+          }
+        io.readExpecting(')');
+        if(isPositions)
+          {
+          // skip label celli
+          io.readLabel();
+          }
+        io.readExpecting('}');
+        return;
+        }
+
+      if(currToken != '(')
+        {
+        throw vtkFoamError() << "Expected (, found " << currToken;
+        }
+      if(io.format() == vtkFoamIOobject::ASCII)
+        {
+        if(size > 0) // avoid invalid access to GetPointer(0)
+          {
+          double* vectorComponentI = vectorListPtr_->GetPointer(0);
+          for(int i = 0; i < size; i++)
+            {
+            io.readExpecting('(');
+            *vectorComponentI++ = io.readNumber();
+            *vectorComponentI++ = io.readNumber();
+            *vectorComponentI++ = io.readNumber();
+            io.readExpecting(')');
+            if(isPositions)
+              {
+              // skip label celli
+              io.readLabel();
+              }
+            }
+          }
+        }
+      else // binary
+        {
+        if(isPositions) // lagrangian/positions (class Cloud)
+          {
+          if(size > 0) // avoid invalid access to GetPointer()
+            {
+            for(int i = 0, index = 0; i < size; i++, index += 3)
+              {
+              io.readExpecting('(');
+              io.read(reinterpret_cast<char*>(vectorListPtr_
+                ->GetPointer(index)), sizeof(double) * 3);
+
+              // skip label celli, label facei and scalar stepFraction
+              const int dummySize = 2 * sizeof(int) + sizeof(double);
+              char dummyBuf[dummySize];
+              io.read(dummyBuf, dummySize);
+              io.readExpecting(')');
+              }
+            }
+          }
+        else // regular vectorField
+          {
+          if(size > 0) // avoid invalid access to GetPointer()
+            {
+            io.read(reinterpret_cast<char*>(vectorListPtr_->GetPointer(0)),
+              size * sizeof(double) * 3);
+            }
+          }
+        }
+      io.readExpecting(')');
+      }
+    else if(currToken == '(')
+      {
+      while(io.read(currToken) && currToken != ')')
+        {
+        if(currToken != '(')
+          {
+          throw vtkFoamError() << "Expected (, found " << currToken;
+          }
+        double v[3];
+        v[0] = io.readNumber();
+        v[1] = io.readNumber();
+        v[2] = io.readNumber();
+        vectorListPtr_->InsertNextTuple(v);
+        io.readExpecting(')');
+        }
+      vectorListPtr_->Squeeze();
+      }
+    else
+      {
+      throw vtkFoamError() << "Expected integer or (, found " << currToken;
+      }
+  }
+
+  const vtkFoamEntryValue& operator>>(vtkStdString& str) const
+  {
+    if(type_ == STRING)
+      {
+      str = vtkStdString(t_.stringToken());
+      }
+    else
+      {
+      //vtkErrorMacro("token type does not match");
+      str = "";
+      }
+    return *this;
+  }
+
+  const vtkFoamEntryValue& operator>>(double& d) const
+  {
+    if(type_ == SCALAR || type_ == LABEL)
+      {
+      d = t_.number();
+      }
+    else
+      {
+      //vtkErrorMacro("token type does not match");
+      d = 0.0;
+      }
+    return *this;
+  }
+
+  const vtkFoamEntryValue& operator>>(int& i) const
+  {
+    if(type_ == LABEL)
+      {
+      i = t_.labelToken();
+      }
+    else
+      {
+      //vtkErrorMacro("token type does not match");
+      i = 0;
+      }
+    return *this;
+  }
+};
+
+// a class that represents an entry of a dictionary. note that an
+// entry can have more than one value.
+struct vtkOpenFOAMReader::vtkFoamEntry
+{
+private:
+  vtkStdString keyword_;
+  vtkstd::vector<vtkFoamEntryValue*> valuePtrs_;
+
+public:
+  vtkFoamEntry() {}
+  ~vtkFoamEntry()
+  {
+    for(size_t i = 0; i < valuePtrs_.size(); i++)
+      {
+      delete valuePtrs_[i];
+      }
+  }
+
+  vtkStdString& keyword() { return keyword_; }
+  size_t size() const { return valuePtrs_.size(); }
+  // returns false if the number of the values is 0 to simplify things
+  bool found() const { return valuePtrs_.size() > 0; }
+  vtkFoamEntryValue& firstValue() const { return *valuePtrs_[0]; }
+  const vtkstd::vector<int>& labelList() const
+  { return firstValue().labelList(); }
+  const vtkstd::vector<vtkstd::vector<int> >& labelListList() const
+  { return firstValue().labelListList(); }
+  vtkDoubleArray& scalarList()
+  { return firstValue().scalarList(); }
+  vtkDoubleArray& vectorList()
+  { return firstValue().vectorList(); }
+  vtkFoamDict& dictionary() // not using firstValue() for breaking constness
+  { return valuePtrs_[0]->dictionary(); }
+  void *ptr() { return firstValue().ptr(); }
+
+  void readDictionary(vtkFoamIOobject& io)
+  {
+    valuePtrs_.push_back(new vtkFoamEntryValue);
+    valuePtrs_.back()->readDictionary(io);
+  }
+
+  // read values of an entry
+  void read(vtkFoamIOobject& io);
+
+  const vtkFoamEntry& operator>>(vtkStdString& str) const
+  {
+    if(!found() || valuePtrs_[0]->type() != vtkFoamEntryValue::STRING)
+      {
+      str = "";
+      }
+    else
+      {
+      valuePtrs_[0]->operator>>(str);
+      }
+    return *this;
+  }
+
+  const vtkFoamEntry& operator>>(double& d) const
+  {
+    if(!found() || (valuePtrs_[0]->type() != vtkFoamEntryValue::SCALAR
+      && valuePtrs_[0]->type() != vtkFoamEntryValue::LABEL))
+      {
+      d = 0.0;
+      }
+    else
+      {
+      valuePtrs_[0]->operator>>(d);
+      }
+    return *this;
+  }
+
+  const vtkFoamEntry& operator>>(int& i) const
+  {
+    if(!found() || valuePtrs_[0]->type() != vtkFoamEntryValue::LABEL)
+      {
+      i = 0;
+      }
+    else
+      {
+      valuePtrs_[0]->operator>>(i);
+      }
+    return *this;
+  }
+};
+
+// a class that holds a FoamFile data structure
+struct vtkOpenFOAMReader::vtkFoamDict: public vtkFoamEntryValue
+{
+public:
+  enum dictType
+    { UNDEFINED, DICTIONARY, SCALARLIST, VECTORLIST, LABELLIST, LABELLISTLIST,
+      UNIFORMLABELLIST, UNIFORMSCALARLIST };
+
+private:
+  dictType type_;
+  vtkstd::vector<vtkFoamEntry*> entryPtrs_;
+  vtkFoamEntry* dummyEntryPtr_;
+
+  vtkFoamDict(const vtkFoamDict &);
+
+public:
+  vtkFoamDict(): vtkFoamEntryValue(), type_(UNDEFINED), dummyEntryPtr_(NULL) {}
+  ~vtkFoamDict()
+  {
+    if(type_ == DICTIONARY)
+      {
+      for(size_t i = 0; i < entryPtrs_.size(); i++)
+        {
+        delete entryPtrs_[i];
+        }
+      }
+    if(dummyEntryPtr_ != NULL)
+      {
+      delete dummyEntryPtr_;
+      }
+  }
+
+  dictType type() const { return type_; }
+  size_t size() const { return entryPtrs_.size(); }
+  vtkFoamEntry& entry(const int i) { return *entryPtrs_[i]; }
+  vtkFoamEntry& lookup(const vtkStdString& keyword)
+  {
+    if(type_ == DICTIONARY)
+      {
+      for(size_t i = 0; i < entryPtrs_.size(); i++)
+        {
+        if(entryPtrs_[i]->keyword() == keyword) // found
+          {
+          return *entryPtrs_[i];
+          }
+        }
+      }
+
+    // not found
+    if(dummyEntryPtr_ == NULL)
+      {
+      dummyEntryPtr_ = new vtkFoamEntry;
+      }
+    return *dummyEntryPtr_;
+  }
+
+  // reads a FoamFile or a subdictionary. if the stream to be read is
+  // a subdictionary the preceding '{' is assumed to have already been
+  // thrown away.
+  bool read(vtkFoamIOobject& io, const bool isSubDictionary = false,
+    const vtkStdString& firstKeyword = "")
+  {
+    try
+      {
+      vtkFoamToken currToken;
+      if(firstKeyword == "")
+        {
+        if(!isSubDictionary)
+          {
+          // polyMesh/points, lagrangian vectors
+          if(io.className() == "vectorField")
+            {
+            readVectorList(io, false);
+            type_ = VECTORLIST;
+            return true;
+            }
+          else if(io.className() == "scalarField") // lagrangian scalars
+            {
+            readScalarList(io);
+            type_ = SCALARLIST;
+            return true;
+            }
+          else if(io.className() == "Cloud") // lagrangian/positions
+            {
+            readVectorList(io, true);
+            type_ = VECTORLIST;
+            return true;
+            }
+          else if(io.className() == "faceList") // polyMesh/faces
+            {
+            readLabelListList(io);
+            type_ = LABELLISTLIST;
+            return true;
+            }
+          else if(io.className() == "labelList") // polyMesh/{owner|neighbour}
+            {
+            readLabelList(io);
+            type_ = LABELLIST;
+            return true;
+            }
+          }
+
+        // read the first token
+        if(!io.read(currToken))
+          {
+          throw vtkFoamError() << "Unexpected EOF";
+          }
+
+        // list of dictionaries is read as a usual dictionary
+        // polyMesh/boundary, point/face/cell-Zones
+        if(!isSubDictionary && currToken.type() == vtkFoamToken::LABEL)
+          {
+          io.readExpecting('(');
+          if(currToken.labelToken() > 0)
+            {
+            if(!io.read(currToken))
+              {
+              throw vtkFoamError() << "Unexpected EOF";
+              }
+            // continue to read as a usual dictionary
+            }
+          else // return as empty dictionary
+            {
+            io.readExpecting(')');
+            type_ = DICTIONARY;
+            return true;
+            }
+          }
+        // some boundary files does not have the number of boundary
+        // patches (e.g. settlingFoam/tank3D). in this case we need to
+        // explicitly read the file as a dictionary.
+        else if(!isSubDictionary && currToken == '('
+          && io.className() == "polyBoundaryMesh") // polyMesh/boundary
+          {
+	    if(!io.read(currToken)) // read the first keyword
+            {
+            throw vtkFoamError() << "Unexpected EOF";
+            }
+          if(currToken == ')') // return as empty dictionary
+            {
+            type_ = DICTIONARY;
+            return true;
+            }
+          }
+	// the following two else-if clauses are for an exceptional
+	// expression of `LABEL{LABELorSCALAR}' without type prefix
+	// (e. g. `2{-0}' in mixedRhoE B.C. in
+	// rhopSonicFoam/shockTube)
+	else if(isSubDictionary && currToken.type() == vtkFoamToken::LABEL)
+	  {
+          token() = currToken.labelToken();
+          type_ = UNIFORMLABELLIST;
+          io.readExpecting('}');
+          return true;
+          }
+	else if(isSubDictionary && currToken.type() == vtkFoamToken::SCALAR)
+	  {
+          token() = currToken.number();
+          type_ = UNIFORMSCALARLIST;
+          io.readExpecting('}');
+          return true;
+          }
+        // return as empty dictionary
+        else if(isSubDictionary && currToken == '}')
+          {
+          type_ = DICTIONARY;
+          return true;
+          }
+        }
+      // if firstKeyword is set read the following stream as subdictionary
+      else
+        {
+        entryPtrs_.push_back(new vtkFoamEntry);
+        entryPtrs_.back()->keyword() = firstKeyword;
+        entryPtrs_.back()->readDictionary(io);
+        if(!io.read(currToken) || currToken == '}' || currToken == ')')
+          {
+          type_ = DICTIONARY;
+          return true;
+          }
+        }
+
+      if(currToken.type() == vtkFoamToken::STRING) // general dictionary
+        {
+        // based on src/OpenFOAM/db/dictionary/dictionaryIO.C
+        do
+          {
+          if(currToken != ';') // ignore empty entry
+            {
+            entryPtrs_.push_back(new vtkFoamEntry);
+            entryPtrs_.back()->keyword() = currToken.stringToken();
+            entryPtrs_.back()->read(io);
+            if(currToken == "FoamFile")
+              {
+              // delete the FoamFile header subdictionary entry
+              delete entryPtrs_.back();
+              entryPtrs_.pop_back();
+              }
+            }
+          } while(io.read(currToken)
+            && (currToken.type() == vtkFoamToken::STRING || currToken == ';'));
+
+        if(currToken.type() == vtkFoamToken::ERROR || currToken == '}'
+          || currToken == ')')
+          {
+          type_ = DICTIONARY;
+          return true;
+          }
+        throw vtkFoamError()
+          << "Expected keyword, closing brace, ; or EOF, found " << currToken;
+        }
+      throw vtkFoamError() << "Bad first keyword for dictionary " << currToken;
+      }
+    catch(vtkFoamError& e)
+      {
+      if(isSubDictionary)
+        {
+        throw;
+        }
+      else
+        {
+        io.setError(e);
+        return false;
+        }
+      }
+  }
+};
+
+void vtkOpenFOAMReader::vtkFoamIOobject::readHeader()
+{
+  vtkFoamToken firstToken;
+
+  readExpecting("FoamFile");
+  readExpecting('{');
+
+  vtkFoamDict headerDict;
+  headerDict.read(*this, true); // throw exception in case of error
+
+  vtkFoamEntry& formatEntry = headerDict.lookup("format");
+  if(!formatEntry.found())
+    {
+    throw vtkFoamError()
+      << "the format entry (binary/ascii) not found in FoamFile header";
+    }
+  vtkStdString format;
+  formatEntry >> format;
+  // case does matter (e. g. "BINARY" is treated as ascii)
+  // see src/OpenFOAM/db/IOstreams/IOstreams/IOstream.C
+  if(format == "binary")
+    {
+    format_ = BINARY;
+    }
+  else
+    {
+    format_ = ASCII;
+    }
+
+  vtkFoamEntry& classEntry = headerDict.lookup("class");
+  if(!classEntry.found())
+    {
+    throw vtkFoamError() << "class name not found in FoamFile header";
+    }
+  classEntry >> headerClassName_;
+
+  vtkFoamEntry& objectEntry = headerDict.lookup("object");
+  if(!objectEntry.found())
+    {
+    throw vtkFoamError() << "object name not found in FoamFile header";
+    }
+  objectEntry >> objectName_;
+}
+
+void vtkOpenFOAMReader::vtkFoamEntryValue::clear()
+{
+  if(managed_)
+    {
+    if(type_ == LABELLIST)
+      {
+      delete labelListPtr_;
+      }
+    else if(type_ == LABELLISTLIST)
+      {
+      delete labelListListPtr_;
+      }
+    else if(type_ == SCALARLIST)
+      {
+      scalarListPtr_->Delete();
+      }
+    else if(type_ == VECTORLIST)
+      {
+      vectorListPtr_->Delete();
+      }
+    else if(type_ == STRINGLIST)
+      {
+      delete stringListPtr_;
+      }
+    else if(type_ == ENTRYVALUELIST)
+      {
+      for(size_t i = 0; i < entryValuePtrs_.size() ; i++)
+        {
+        delete entryValuePtrs_[i];
+        }
+      }
+    else if(type_ == DICTIONARY)
+      {
+      delete dictPtr_;
+      }
+    }
+}
+
+// general-purpose list reader - guess the type of the list and read
+// it. only supports ascii format and assumes the preceding '(' has
+// already been thrown away.  the reader supports nested list with
+// variable lengths (e. g. `((token token) (token token token)).'
+// also supports compound of tokens and lists (e. g. `((token token)
+// token)') only if a list comes as the first value.
+void vtkOpenFOAMReader::vtkFoamEntryValue::readList(vtkFoamIOobject& io)
+{
+  io.read(t_);
+
+  // initial guess of the list type
+  if(t_.type() == vtkFoamToken::LABEL)
+    {
+    // if the first token is of type LABEL it might be either an element of
+    // a labelList or the size of a sublist so proceed to the next token
+    vtkFoamToken nextToken;
+    if(!io.read(nextToken))
+      {
+      throw vtkFoamError() << "Unexpected EOF";
+      }
+    if(nextToken.type() == vtkFoamToken::LABEL)
+      {
+      labelListPtr_ = new intVector;
+      labelListPtr_->value.push_back(t_.labelToken());
+      labelListPtr_->value.push_back(nextToken.labelToken());
+      type_ = LABELLIST;
+      }
+    else if(nextToken.type() == vtkFoamToken::SCALAR)
+      {
+      scalarListPtr_ = vtkDoubleArray::New();
+      scalarListPtr_->InsertNextValue(t_.number());
+      scalarListPtr_->InsertNextValue(nextToken.number());
+      type_ = SCALARLIST;
+      }
+    else if(nextToken == '(') // list of list: read recursively
+      {
+      entryValuePtrs_.push_back(new vtkFoamEntryValue);
+      entryValuePtrs_.back()->readList(io);
+      type_ = ENTRYVALUELIST;
+      }
+    else if(nextToken == ')') // list with only one label element
+      {
+      labelListPtr_ = new intVector;
+      labelListPtr_->value.push_back(t_.labelToken());
+      type_ = LABELLIST;
+      return;
+      }
+    else
+      {
+      throw vtkFoamError() << "Expected number, ( or ), found " << nextToken;
+      }
+    }
+  else if(t_.type() == vtkFoamToken::SCALAR)
+    {
+    scalarListPtr_ = vtkDoubleArray::New();
+    scalarListPtr_->InsertNextValue(t_.number());
+    type_ = SCALARLIST;
+    }
+  // if the first word is a string we have to read another token to determine
+  // if the first word is a keyword for the following dictionary
+  else if(t_.type() == vtkFoamToken::STRING)
+    {
+    vtkFoamToken nextToken;
+    if(!io.read(nextToken))
+      {
+      throw vtkFoamError() << "Unexpected EOF";
+      }
+    if(nextToken.type() == vtkFoamToken::STRING) // list of strings
+      {
+      stringListPtr_ = new stringVector;
+      stringListPtr_->value.push_back(t_.stringToken());
+      stringListPtr_->value.push_back(nextToken.stringToken());
+      type_ = STRINGLIST;
+      }
+    // dictionary with the already read stringToken as the first keyword
+    else if(nextToken == '{')
+      {
+      if(t_.stringToken() == "")
+        {
+        throw "Empty string is invalid as a keyword for dictionary entry";
+        }
+      readDictionary(io, t_.stringToken());
+      // the dictionary read as list has the entry terminator ';' so
+      // we have to skip it
+      return;
+      }
+    else if(nextToken == ')') // list with only one string element
+      {
+      stringListPtr_ = new stringVector;
+      stringListPtr_->value.push_back(t_.stringToken());
+      type_ = STRINGLIST;
+      return;
+      }
+    else
+      {
+      throw vtkFoamError() << "Expected string, { or ), found " << nextToken;
+      }
+    }
+  else if(t_ == '(') // list of lists: read recursively
+    {
+    entryValuePtrs_.push_back(new vtkFoamEntryValue);
+    entryValuePtrs_.back()->readList(io);
+    // read all the following values as arbitrary entryValues
+    // the alphaContactAngle b.c. in multiphaseInterFoam/damBreak4phase
+    // reaquires this treatment (reading by readList() is not enough)
+    do
+      {
+      entryValuePtrs_.push_back(new vtkFoamEntryValue);
+      entryValuePtrs_.back()->read(io);
+      }
+    while(*entryValuePtrs_.back() != ')' && *entryValuePtrs_.back() != '}'
+      && *entryValuePtrs_.back() != ';');
+
+    if(*entryValuePtrs_.back() != ')')
+      {
+      throw vtkFoamError() << "Expected ) before "
+        << entryValuePtrs_.back()->token();
+      }
+
+    // delete ')'
+    delete entryValuePtrs_.back();
+    entryValuePtrs_.pop_back();
+    type_ = ENTRYVALUELIST;
+    return;
+    }
+  else if(t_ == ')') // empty list
+    {
+    type_ = EMPTYLIST;
+    return;
+    }
+
+  while(io.read(t_) && t_ != ')')
+    {
+    if(type_ == LABELLIST)
+      {
+      if(t_.type() == vtkFoamToken::SCALAR) // switch to scalarList
+        {
+        // labelListPtr_ and scalarListPtr_ are packed into a single union so
+        // we need a temprary pointer
+        vtkDoubleArray* slPtr = vtkDoubleArray::New();
+        const vtkIdType size = labelListPtr_->value.size();
+        slPtr->SetNumberOfValues(size);
+        for(vtkIdType i = 0; i < size; i++)
+          {
+          slPtr->SetValue(i, double(labelListPtr_->value[i]));
+          }
+        delete labelListPtr_;
+        slPtr->InsertNextValue(t_.number());
+        scalarListPtr_ = slPtr; // copy after labelListPtr_ is deleted
+        type_ = SCALARLIST;
+        }
+      else if(t_.type() == vtkFoamToken::LABEL)
+        {
+        labelListPtr_->value.push_back(t_.labelToken());
+        }
+      else
+        {
+        throw vtkFoamError() << "Expected a number, found " << t_;
+        }
+      }
+    else if(type_ == SCALARLIST)
+      {
+      if(t_.isNumber())
+        {
+        scalarListPtr_->InsertNextValue(t_.number());
+        }
+      else
+        {
+        throw vtkFoamError() << "Expected a number, found " << t_;
+        }
+      }
+    else if(type_ == STRINGLIST)
+      {
+      if(t_.type() == vtkFoamToken::STRING)
+        {
+        stringListPtr_->value.push_back(t_.stringToken());
+        }
+      else
+        {
+        throw vtkFoamError() << "Expected a string, found " << t_;
+        }
+      }
+    else if(type_ == ENTRYVALUELIST)
+      {
+      if(t_.type() == vtkFoamToken::LABEL)
+        {
+        // skip the number of elements to make things simple
+        if(!io.read(t_))
+          {
+          throw vtkFoamError() << "Unexpected EOF";
+          }
+        }
+      if(t_ != '(')
+        {
+        throw vtkFoamError() << "Expected (, found " << t_;
+        }
+      entryValuePtrs_.push_back(new vtkFoamEntryValue);
+      entryValuePtrs_.back()->readList(io);
+      }
+    else
+      {
+      throw vtkFoamError() << "Unexpected token " << t_;
+      }
+    }
+
+  if(type_ == SCALARLIST)
+    {
+    scalarListPtr_->Squeeze();
+    }
+}
+
+// a list of dictionaries is actually read as a dictionary
+void vtkOpenFOAMReader::vtkFoamEntryValue::readDictionary(vtkFoamIOobject& io,
+  const vtkStdString& firstKeyword)
+{
+  dictPtr_ = new vtkFoamDict;
+  type_ = DICTIONARY;
+  dictPtr_->read(io, true, firstKeyword);
+}
+
+// guess the type of the given entry value and read it
+void vtkOpenFOAMReader::vtkFoamEntryValue::read(vtkFoamIOobject& io)
+{
+  if(!io.read(t_))
+    {
+    throw vtkFoamError() << "Unexpected EOF";
+    }
+
+  if(t_ == '{')
+    {
+    readDictionary(io);
+    return;
+    }
+  // for reading sublist from vtkFoamEntryValue::readList() or there
+  // are cases where lists without the (non)uniform keyword appear
+  // (e. g. coodles/pitsDaily/0/U, Hrv's uniformFixedValue b.c.)
+  else if(t_ == '(')
+    {
+    readList(io);
+    return;
+    }
+  else if(t_ == "uniform")
+    {
+    if(!io.read(t_))
+      {
+      throw vtkFoamError()
+        << "Expected a uniform value or a list, found unexpected EOF";
+      }
+    if(t_ == '(')
+      {
+      readList(io);
+      }
+    else if(t_.type() == vtkFoamToken::LABEL)
+      {
+      type_ = LABEL;
+      }
+    else if(t_.type() == vtkFoamToken::SCALAR)
+      {
+      type_ = SCALAR;
+      }
+    else if(t_.type() == vtkFoamToken::STRING)
+      {
+      type_ = STRING;
+      }
+    else // unexpected punctuation token
+      {
+      throw vtkFoamError() << "Expected number, string or (, found " << t_;
+      }
+    isUniform_ = true;
+    }
+  else if(t_ == "nonuniform")
+    {
+    if(!io.read(t_))
+      {
+      throw vtkFoamError() << "Expected list type specifier, found EOF";
+      }
+    if(t_ == "List<vector>")
+      {
+      isUniform_ = false;
+      readVectorList(io, false);
+      }
+    else if(t_ == "List<scalar>")
+      {
+      isUniform_ = false;
+      readScalarList(io);
+      }
+    // List<bool> is read as List<label>
+    else if(t_ == "List<label>" || t_ == "List<bool>")
+      {
+      isUniform_ = false;
+      readLabelList(io);
+      }
+    // an empty list doesn't have a list type specifier
+    else if(t_.type() == vtkFoamToken::LABEL && t_.labelToken() == 0)
+      {
+      type_ = EMPTYLIST;
+      isUniform_ = false;
+      io.readExpecting('(');
+      io.readExpecting(')');
+      }
+    else
+      {
+      throw vtkFoamError() << "Unsupported nonuniform list type " << t_;
+      }
+    }
+  // zones have list without a uniform/nonuniform keyword
+  // List<bool> is read as List<label>
+  // (e. g. flipMap entry in point/face/cellZones)
+  else if(t_ == "List<label>" || t_ == "List<bool>")
+    {
+    isUniform_ = false;
+    readLabelList(io);
+    }
+  else if(t_.type() == vtkFoamToken::LABEL)
+    {
+    type_ = LABEL;
+    }
+  else if(t_.type() == vtkFoamToken::SCALAR)
+    {
+    type_ = SCALAR;
+    }
+  else if(t_.type() == vtkFoamToken::STRING)
+    {
+    type_ = STRING;
+    }
+  else if(t_.type() == vtkFoamToken::PUNCTUATION)
+    {
+    type_ = PUNCTUATION;
+    }
+}
+
+// read values of an entry
+void vtkOpenFOAMReader::vtkFoamEntry::read(vtkFoamIOobject& io)
 {
-  vtkstd::vector< int > value;
-};
+  for(;;)
+    {
+    valuePtrs_.push_back(new vtkFoamEntryValue);
+    valuePtrs_.back()->read(io);
 
-struct intVectorVector
-{
-  vtkstd::vector< vtkstd::vector< int > > value;
-};
+    if(valuePtrs_.size() >= 2)
+      {
+      vtkFoamEntryValue& secondLastValue
+        = *valuePtrs_[valuePtrs_.size() - 2];
+      if(secondLastValue.type() == vtkFoamEntryValue::LABEL)
+        {
+        vtkFoamEntryValue& lastValue = *valuePtrs_.back();
+
+        // a zero-sized nonuniform list without prefixing "nonuniform"
+        // keyword nor list type specifier (i. e. `0()';
+        // e. g. simpleEngine/0/polyMesh/pointZones) requires special
+        // care (one with nonuniform prefix is treated within
+        // vtkFoamEntryValue::read()). still this causes errornous
+        // behavior for `0 nonuniform 0()' but this should be extremely
+        // rare
+        if(lastValue.type() == vtkFoamEntryValue::EMPTYLIST
+          && secondLastValue == 0)
+          {
+          delete valuePtrs_.back();
+          valuePtrs_.pop_back(); // delete the last value
+          valuePtrs_.back()->setEmptyList(); // mark new last value as empty
+          }
+        // for an exceptional expression of `LABEL{LABELorSCALAR}' without
+        // type prefix (e. g. `2{-0}' in mixedRhoE B.C. in
+        // rhopSonicFoam/shockTube)
+        else if(lastValue.type() == vtkFoamEntryValue::DICTIONARY)
+          {
+          if(lastValue.dictionary().type() == vtkFoamDict::UNIFORMLABELLIST)
+            {
+            const int size = secondLastValue.token().labelToken();
+            const int value = lastValue.dictionary().token().labelToken();
+            // delete last two values
+            delete valuePtrs_.back();
+            valuePtrs_.pop_back();
+            delete valuePtrs_.back();
+            valuePtrs_.pop_back();
+            // make new labelList
+            valuePtrs_.push_back(new vtkFoamEntryValue);
+            valuePtrs_.back()->makeLabelList(value, size);
+            }
+          else if(lastValue.dictionary().type()
+            == vtkFoamDict::UNIFORMSCALARLIST)
+            {
+            const int size = secondLastValue.token().labelToken();
+            const double value = lastValue.dictionary().token().number();
+            // delete last two values
+            delete valuePtrs_.back();
+            valuePtrs_.pop_back();
+            delete valuePtrs_.back();
+            valuePtrs_.pop_back();
+            // make new labelList
+            valuePtrs_.push_back(new vtkFoamEntryValue);
+            valuePtrs_.back()->makeScalarList(value, size);
+            }
+          }
+        }
+      }
 
-struct faceVectorVector
-{
-  vtkstd::vector< vtkstd::vector< face > > value;
-};
+    if(*valuePtrs_.back() == ';')
+      {
+      delete valuePtrs_.back();
+      valuePtrs_.pop_back();
+      break;
+      }
+    else if(valuePtrs_.back()->type() == vtkFoamEntryValue::DICTIONARY)
+      {
+      // subdictionary is not suffixed by an entry terminator ';'
+      break;
+      }
+    else if(*valuePtrs_.back() == '}' || *valuePtrs_.back() == ')')
+      {
+      throw vtkFoamError() << "Unmatched " << valuePtrs_.back()->token();
+      }
+    }
+}
 
 vtkOpenFOAMReader::vtkOpenFOAMReader()
 {
@@ -95,60 +2168,144 @@
   this->FileName = NULL;
 
   //VTK CLASSES
-  this->Points = vtkPoints::New();
+  this->PatchDataArraySelection = vtkDataArraySelection::New();
   this->CellDataArraySelection = vtkDataArraySelection::New();
+  this->PointDataArraySelection = vtkDataArraySelection::New();
+
+  // Setup the Selection observer for the above selection arrays
+  this->SelectionObserver = vtkCallbackCommand::New();
+  this->SelectionObserver
+    ->SetCallback(&vtkOpenFOAMReader::SelectionModifiedCallback);
+  this->SelectionObserver->SetClientData(this);
+
+  this->PatchDataArraySelection
+    ->AddObserver(vtkCommand::ModifiedEvent,this->SelectionObserver);
+  this->CellDataArraySelection
+    ->AddObserver(vtkCommand::ModifiedEvent,this->SelectionObserver);
+  this->PointDataArraySelection
+    ->AddObserver(vtkCommand::ModifiedEvent,this->SelectionObserver);
+
+  // Initialise the Selection status arrays
+  this->CellSelectionOldStatus = 0;
+  this->PointSelectionOldStatus = 0;
+  this->PatchSelectionOldStatus = 0;
+  this->CellSelectionStatus = 0;
+  this->PointSelectionStatus = 0;
+  this->PatchSelectionStatus = 0;
 
   //DATA COUNTS
-  this->NumFaces = 0;
-  this->NumPoints = 0;
   this->NumCells = 0;
 
   this->TimeStepData = new stringVector;
-  this->Path = new stdString;
-  this->PathPrefix = new stdString;
+  this->LagrangianTimeStepData = new stringVector;
+  this->OldFileName = new vtkStdString;
+  this->PathPrefix = new vtkStdString;
   this->PolyMeshPointsDir = new stringVector;
   this->PolyMeshFacesDir = new stringVector;
-  this->BoundaryNames = new stringVector;
-  this->PointZoneNames = new stringVector;
-  this->FaceZoneNames = new stringVector;
-  this->CellZoneNames = new stringVector;
-
-  this->FacePoints = new intVectorVector;
-  this->FacesOwnerCell = new intVectorVector;
-  this->FacesNeighborCell = new intVectorVector;
-  this->FacesOfCell = new faceVectorVector;
-  this->SizeOfBoundary = new intVector;
 
   //DATA TIMES
   this->NumberOfTimeSteps = 0;
   this->Steps = NULL;
+  this->TimeNames = vtkStringArray::New();
   this->TimeStep = 0;
   this->TimeStepRange[0] = 0;
   this->TimeStepRange[1] = 0;
-  this->RequestInformationFlag = true;
+
+  // for caching mesh
+  this->CacheMesh = 1;
+  this->OldTimeStep = -1;
+  this->InternalMesh = NULL;
+  this->BoundaryMesh = NULL;
+  this->BoundaryDict = NULL;
+  this->FaceOwner = NULL;
+  this->PointZoneMesh = NULL;
+  this->FaceZoneMesh = NULL;
+  this->CellZoneMesh = NULL;
+
+  // for reading zones
+  this->ReadZones = 0; // turned off by default
+  // for accumulating patches across time-steps (Turned on by default)
+  // As of now, only accumulation of patches is implemented
+  this->KeepPatches = 1;
 }
 
 vtkOpenFOAMReader::~vtkOpenFOAMReader()
 {
   vtkDebugMacro(<<"DeConstructor");
-  this->Points->Delete();
+
+  // Delete the Observers before deleting the Selection Arrays !!!
+  this->PatchDataArraySelection->RemoveObserver(this->SelectionObserver);
+  this->CellDataArraySelection->RemoveObserver(this->SelectionObserver);
+  this->PointDataArraySelection->RemoveObserver(this->SelectionObserver);
+
+  // Now delete the Selection Observer itself
+  this->SelectionObserver->Delete();
+
+  // Finally delete the Selection Arrays
+  this->PatchDataArraySelection->Delete();
   this->CellDataArraySelection->Delete();
-  delete [] this->Steps;
+  this->PointDataArraySelection->Delete();
 
-  delete this->TimeStepData;
-  delete this->Path;
+  delete this->OldFileName;
   delete this->PathPrefix;
   delete this->PolyMeshPointsDir;
   delete this->PolyMeshFacesDir;
-  delete this->BoundaryNames;
-  delete this->PointZoneNames;
-  delete this->FaceZoneNames;
-  delete this->CellZoneNames;
-  delete this->FacePoints;
-  delete this->FacesOwnerCell;
-  delete this->FacesNeighborCell;
-  delete this->FacesOfCell;
-  delete this->SizeOfBoundary;
+  this->TimeNames->Delete();
+  delete this->TimeStepData;
+  delete this->LagrangianTimeStepData;
+
+  this->ClearMeshes();
+}
+
+void vtkOpenFOAMReader::ClearMeshes()
+{
+  if(this->FaceOwner != NULL)
+    {
+    this->FaceOwner->Delete();
+    this->FaceOwner = NULL;
+    }
+  if(this->InternalMesh != NULL)
+    {
+    this->InternalMesh->Delete();
+    this->InternalMesh = NULL;
+    }
+  if(this->BoundaryMesh != NULL)
+    {
+    for(size_t i = 0; i < this->BoundaryMesh->value.size(); i++)
+      {
+      this->BoundaryMesh->value[i]->Delete();
+      }
+    delete this->BoundaryMesh;
+    this->BoundaryMesh = NULL;
+    }
+  delete this->BoundaryDict;
+  if(this->PointZoneMesh != NULL)
+    {
+    for(size_t i = 0; i < this->PointZoneMesh->value.size(); i++)
+      {
+      this->PointZoneMesh->value[i]->Delete();
+      }
+    delete this->PointZoneMesh;
+    this->PointZoneMesh = NULL;
+    }
+  if(this->FaceZoneMesh != NULL)
+    {
+    for(size_t i = 0; i < this->FaceZoneMesh->value.size(); i++)
+      {
+      this->FaceZoneMesh->value[i]->Delete();
+      }
+    delete this->FaceZoneMesh;
+    this->FaceZoneMesh = NULL;
+    }
+  if(this->CellZoneMesh != NULL)
+    {
+    for(size_t i = 0; i < this->CellZoneMesh->value.size(); i++)
+      {
+      this->CellZoneMesh->value[i]->Delete();
+      }
+    delete this->CellZoneMesh;
+    this->CellZoneMesh = NULL;
+    }
 }
 
 int vtkOpenFOAMReader::RequestData(
@@ -160,13 +2317,55 @@
   vtkInformation* outInfo = outputVector->GetInformationObject(0);
   vtkMultiBlockDataSet *output = vtkMultiBlockDataSet::SafeDownCast(
     outInfo->Get(vtkMultiBlockDataSet::DATA_OBJECT()));
+
   if(!this->FileName)
     {
     vtkErrorMacro("FileName has to be specified!");
     return 0;
     }
-  this->CreateDataSet(output);
-  return 1;
+
+  if(this->NumberOfTimeSteps == 0)
+    {
+    vtkErrorMacro(<< this->FileName << " contains no timestep data.");
+    return 0;
+    }
+    
+  if(this->TimeStep < this->TimeStepRange[0]
+     || this->TimeStep > this->TimeStepRange[1])
+    {
+    vtkErrorMacro("TimeStep out of range");
+    return 0;
+    }
+
+#if PARAVIEW_VERSION_MAJOR >= 3
+  if(outInfo->Has(vtkStreamingDemandDrivenPipeline::UPDATE_TIME_STEPS()))
+    {
+    double* requestedTimeValues
+      = outInfo->Get(vtkStreamingDemandDrivenPipeline::UPDATE_TIME_STEPS());
+    int nSteps
+      = outInfo->Length(vtkStreamingDemandDrivenPipeline::TIME_STEPS());
+    double* steps
+      = outInfo->Get(vtkStreamingDemandDrivenPipeline::TIME_STEPS());
+    int timeI = 0;
+    while(timeI < nSteps - 1 && steps[timeI] < requestedTimeValues[0])
+      {
+      timeI++;
+      }
+    this->SetTimeStep(timeI);
+
+    outInfo->Set(vtkDataObject::DATA_TIME_STEPS(), &steps[timeI], 1);
+    }
+#endif
+  this->MakeTimeStepData();
+
+  const int ret = this->CreateDataSet(output, this->TimeStep);
+
+  // update selection status
+  this->CellSelectionOldStatus = this->CellSelectionStatus;
+  this->PointSelectionOldStatus = this->PointSelectionStatus;
+  this->PatchSelectionOldStatus = this->PatchSelectionStatus;
+
+  return ret;
 }
 
 void vtkOpenFOAMReader::PrintSelf(ostream& os, vtkIndent indent)
@@ -175,7 +2374,7 @@
   this->Superclass::PrintSelf(os,indent);
   os << indent << "File Name: "
      << (this->FileName ? this->FileName : "(none)") << "\n";
-  os << indent << "Number Of Nodes: " << this->NumPoints << "\n";
+//  os << indent << "Number Of Nodes: " << this->NumPoints << "\n";
   os << indent << "Number Of Cells: " << this->NumCells << "\n";
   os << indent << "Number of Time Steps: " << this->NumberOfTimeSteps << endl;
   os << indent << "TimeStepRange: " 
@@ -190,74 +2389,58 @@
   vtkInformationVector **vtkNotUsed(inputVector),
   vtkInformationVector *outputVector)
 {
-  if(!this->FileName)
+  if(!this->FileName || strlen(this->FileName) == 0)
     {
     vtkErrorMacro("FileName has to be specified!");
     return 0;
     }
-  vtkDebugMacro(<<"Request Info");
-  if(RequestInformationFlag)
+  vtkDebugMacro(<<"Request Info: " <<this->FileName);
+
+  if(*this->OldFileName != vtkStdString(this->FileName))
     {
-    vtkDebugMacro(<<this->FileName);
-    this->Path->value.append(this->FileName);
-    this->ReadControlDict();
+    // updated this->OldFileName
+    *this->OldFileName = vtkStdString(this->FileName);
+
+    // clear prior case information
+    this->OldTimeStep = -1;
+    if(this->Steps != NULL)
+      {
+      delete [] this->Steps;
+      this->Steps = NULL;
+      }
+    this->ClearMeshes();
+
+    this->CellDataArraySelection->RemoveAllArrays();
+    this->PointDataArraySelection->RemoveAllArrays();
+    this->PatchDataArraySelection->RemoveAllArrays();
+
+    // make case information
+    if(!this->ReadControlDict(this->FileName))
+      {
+      return 0;
+      }
+    if(this->NumberOfTimeSteps == 0)
+      {
+      vtkErrorMacro(<< this->FileName << " contains no timestep data.");
+      return 0;
+      }
     this->TimeStepRange[0] = 0;
-    this->TimeStepRange[1] = this->NumberOfTimeSteps-1;
+    this->TimeStepRange[1] = this->NumberOfTimeSteps;
     this->PopulatePolyMeshDirArrays();
     outputVector->GetInformationObject(0)->Set(
       vtkStreamingDemandDrivenPipeline::TIME_STEPS(), this->Steps,
       this->NumberOfTimeSteps);
-    this->RequestInformationFlag = false;
-    }
-
-  //Add scalars and vectors to metadata
-  //create path to current time step
-  vtksys_ios::stringstream tempPath;
-  tempPath << this->PathPrefix->value.c_str();
-  tempPath << this->Steps[this->TimeStep];
-
-  //open the directory and get num of files
-  int numSolvers;
-  vtkDirectory * directory = vtkDirectory::New();
-  int opened = directory->Open(tempPath.str().c_str());
-  if(opened)
-    {
-    numSolvers = directory->GetNumberOfFiles();
-    }
-  else
-    {
-    numSolvers = -1; //no dir
+#if PARAVIEW_VERSION_MAJOR >= 3
+    double timeRange[2];
+    timeRange[0] = this->Steps[0];
+    timeRange[1] = this->Steps[this->NumberOfTimeSteps - 1];
+    outputVector->GetInformationObject(0)->Set(
+      vtkStreamingDemandDrivenPipeline::TIME_RANGE(), timeRange, 2);
+#endif
     }
 
-  //clear prior timestep data
-  this->TimeStepData->value.clear();
-
-  //loop over all files and locate
-  //volScalars and volVectors
-  for(int j = 0; j < numSolvers; j++)
-    {
-    const char * tempSolver = directory->GetFile(j);
-    if(tempSolver != (char *)"polyMesh")
-      {
-      if(tempSolver != (char *)"." && tempSolver != (char *)"..")
-        {
-        vtkstd::string type(this->GetDataType(tempPath.str().c_str(),
-                                              tempSolver));
-        if(strcmp(type.c_str(), "Scalar") == 0)
-        {
-          this->TimeStepData->value.push_back(vtkstd::string(tempSolver));
-          this->CellDataArraySelection->AddArray(tempSolver);
-          }
-        else if(strcmp(type.c_str(), "Vector") == 0)
-          {
-          this->TimeStepData->value.push_back(vtkstd::string(tempSolver));
-          this->CellDataArraySelection->AddArray(tempSolver);
-          }
-        }
-      }
-    }
+  this->MakeTimeStepData();
 
-  directory->Delete();
   return 1;
 }
 
@@ -304,2071 +2487,1860 @@
   return;
 }
 
-// ****************************************************************************
-//  Method: vtkOpenFOAMReader::CombineOwnerNeigbor
-//
-//  Purpose:
-//  add Owner faces to the faces of a cell and then add the neighor faces
-//
-// ****************************************************************************
-void vtkOpenFOAMReader::CombineOwnerNeigbor()
-{
-  vtkDebugMacro(<<"Combine owner & neighbor faces");
-  //reintialize faces of the cells
-  face tempFace;
-  this->FacesOfCell->value.clear();
-  this->FacesOfCell->value.resize(this->NumCells);
+//
+// POINT METHODS
+//
+int vtkOpenFOAMReader::GetNumberOfPointArrays()
+{
+  return this->PointDataArraySelection->GetNumberOfArrays();
+}
+
+const char* vtkOpenFOAMReader::GetPointArrayName(int index)
+{
+  return this->PointDataArraySelection->GetArrayName(index);
+}
+
+int vtkOpenFOAMReader::GetPointArrayStatus(const char* name)
+{
+  return this->PointDataArraySelection->ArrayIsEnabled(name);
+}
+
+void vtkOpenFOAMReader::SetPointArrayStatus(const char* name, int status)
+{
+  if(status)
+    {
+    this->PointDataArraySelection->EnableArray(name);
+    }
+  else
+    {
+    this->PointDataArraySelection->DisableArray(name);
+    }
+  return;
+}
+
+void vtkOpenFOAMReader::DisableAllPointArrays()
+{
+  this->PointDataArraySelection->DisableAllArrays();
+  return;
+}
+
+void vtkOpenFOAMReader::EnableAllPointArrays()
+{
+  this->PointDataArraySelection->EnableAllArrays();
+  return;
+}
+
+
+//
+// PATCH METHODS
+//
+int vtkOpenFOAMReader::GetNumberOfPatchArrays()
+{
+  return this->PatchDataArraySelection->GetNumberOfArrays();
+}
+
+const char* vtkOpenFOAMReader::GetPatchArrayName(int index)
+{
+  return this->PatchDataArraySelection->GetArrayName(index);
+}
+
+int vtkOpenFOAMReader::GetPatchArrayStatus(const char* name)
+{
+  return this->PatchDataArraySelection->ArrayIsEnabled(name);
+}
+
+void vtkOpenFOAMReader::SetPatchArrayStatus(const char* name, int status)
+{
+  if(status)
+    {
+    this->PatchDataArraySelection->EnableArray(name);
+    }
+  else
+    {
+    this->PatchDataArraySelection->DisableArray(name);
+    }
+  return;
+}
+
+void vtkOpenFOAMReader::DisableAllPatchArrays()
+{
+  this->PatchDataArraySelection->DisableAllArrays();
+  return;
+}
+
+void vtkOpenFOAMReader::EnableAllPatchArrays()
+{
+  this->PatchDataArraySelection->EnableAllArrays();
+  return;
+}
+
+// Define the Selection Observer 
+void vtkOpenFOAMReader::SelectionModifiedCallback(vtkObject*,
+  unsigned long, void* clientdata, void*)
+{
+  static_cast<vtkOpenFOAMReader*>(clientdata)->SelectionModified();
+}
+
+void vtkOpenFOAMReader::SelectionModified()
+{
+  // Change the pipeline modification time to force update
+  // Update the selection status to detect changes
+
+  // Cell Selection Arrays
+  // we loop from the last element of the selection array since we'd
+  // like to place newly added variables to extra higher bits (we
+  // wouldn't like to move the locations of already existing variables)
+  this->CellSelectionStatus = 0;
+  for(int i = this->GetNumberOfCellArrays() - 1; i >= 0; i--)
+    {
+    this->CellSelectionStatus = (this->CellSelectionStatus << 1) 
+      + this->GetCellArrayStatus(this->GetCellArrayName(i));
+    }
+  // if the status flag overflows we'd like to fall onto the safe side
+  // (properly refresh the dataset without mesh caching) so we use the
+  // last bit (LSB) as overflow control flag.  if overflow occurs the
+  // LSB is negated from the old status
+  if(this->GetNumberOfCellArrays()
+    <= static_cast<int>(sizeof(this->CellSelectionStatus) * 8 - 1))
+    {
+    this->CellSelectionStatus = (this->CellSelectionStatus << 1)
+      + (this->CellSelectionOldStatus & 0x1);
+    }
+  else
+    {
+    this->CellSelectionStatus = (this->CellSelectionStatus << 1) 
+      + ((~this->CellSelectionOldStatus) & 0x1);
+    }
+  // Point Selection Arrays
+  // this->PointSelectionOldStatus = this->PointSelectionStatus;
+  this->PointSelectionStatus = 0;
+  for(int i = this->GetNumberOfPointArrays() - 1; i >= 0; i--)
+    {
+    this->PointSelectionStatus = (this->PointSelectionStatus << 1) 
+      + this->GetPointArrayStatus(this->GetPointArrayName(i));
+    }
+  if(this->GetNumberOfPointArrays()
+    <= static_cast<int>(sizeof(this->PointSelectionStatus) * 8 - 1))
+    {
+    this->PointSelectionStatus = (this->PointSelectionStatus << 1)
+      + (this->PointSelectionOldStatus & 0x1);
+    }
+  else
+    {
+    this->PointSelectionStatus = (this->PointSelectionStatus << 1) 
+      + (~(this->PointSelectionOldStatus) & 0x1);
+    }
+  // Patch Selection Arrays
+  // this->PatchSelectionOldStatus = this->PatchSelectionStatus;
+  this->PatchSelectionStatus = 0;
+  for(int i = this->GetNumberOfPatchArrays() - 1; i >= 0; i--)
+    {
+    this->PatchSelectionStatus = (this->PatchSelectionStatus << 1) 
+      + this->GetPatchArrayStatus(this->GetPatchArrayName(i));
+    }
+  if(this->GetNumberOfPatchArrays()
+    <= static_cast<int>(sizeof(this->PatchSelectionStatus) * 8 - 1))
+    {
+    this->PatchSelectionStatus = (this->PatchSelectionStatus << 1)
+      + (this->PatchSelectionOldStatus & 0x1);
+    }
+  else
+    {
+    this->PatchSelectionStatus = (this->PatchSelectionStatus << 1) 
+      + (~(this->PatchSelectionOldStatus) & 0x1);
+    }
+
+  // Indicate that the pipeline needs to be updated (VTK Command)
+  this->Modified();
+}
+
+
+// create field data lists and cell/point array selection lists
+void vtkOpenFOAMReader::MakeTimeStepData()
+{
+  // Read the patches from the boundary file into selection array
+  int timeState = this->TimeStep;
+  vtkFoamDict* boundaryDictPtr
+    = this->GatherBlocks("boundary", timeState, true);
+  if(boundaryDictPtr == NULL)
+    {
+    vtkErrorMacro(<< "Couldn't read polyMesh/boundary");
+    return;
+    }
+
+  vtkFoamDict& boundaryDictTmp = *boundaryDictPtr;
+
+  // Add the internal mesh by default always
+  const vtkStdString tmpStr = "Internal Mesh";
+  this->PatchDataArraySelection->AddArray(tmpStr.c_str());
+
+  // iterate through each entry in the boundary file
+  for(size_t i = 0; i < boundaryDictTmp.size(); i++)
+    {
+    vtkFoamEntry& boundaryEntryI = boundaryDictTmp.entry(i);
+    int nFaces;
+    boundaryEntryI.dictionary().lookup("nFaces") >> nFaces;
+
+    // extract name of the current patch for insertion
+    const vtkStdString& boundaryNameI = boundaryEntryI.keyword();
+
+    // If the size of patch becomes zero, but KeepPatches is selected,
+    // do not remove from list, else, remove patch from list
+    if(nFaces == 0)
+      {
+      if((this->PatchDataArraySelection->ArrayExists(boundaryNameI.c_str()))
+        && (!this->GetKeepPatches()))
+        {
+        this->PatchDataArraySelection->RemoveArrayByName(boundaryNameI.c_str());
+        }
+      }
+    else
+      {
+      this->PatchDataArraySelection->AddArray(boundaryNameI.c_str());
+      }
+    }
+  delete boundaryDictPtr;
+
+  // Commented out the "RemoveAllArrays()" line....
+  //this->CellDataArraySelection->RemoveAllArrays();
+
+  //clear prior timestep data
+  this->TimeStepData->value.clear();
+
+  //Add scalars and vectors to metadata
+  //create path to current time step
+  vtkStdString tempPath
+    = *this->PathPrefix + this->TimeNames->GetValue(this->TimeStep);
+
+  //open the directory and get num of files
+  vtkDirectory * directory = vtkDirectory::New();
+  if(!directory->Open(tempPath.c_str()))
+    {
+    // no data
+    directory->Delete();
+    return;
+    }
 
-  //add owner faces to cell
-  for(int i = 0; i < (int)this->FacesOwnerCell->value.size(); i++)
+  //loop over all files and locate
+  //volScalars and volVectors
+  int nFieldFiles = directory->GetNumberOfFiles();
+  for(int j = 0; j < nFieldFiles; j++)
     {
-    for(int j = 0; j < (int)this->FacesOwnerCell->value[i].size(); j++)
+    const vtkStdString fieldFile(directory->GetFile(j));
+    const int len = fieldFile.length();
+
+    // excluded extensions based on src/OpenFOAM/OSspecific/Unix/Unix.C
+    if(!directory->FileIsDirectory(fieldFile.c_str())
+      && fieldFile.substr(len - 1) != "~"
+      && (len < 4 || (fieldFile.substr(len - 4) != ".bak"
+      && fieldFile.substr(len - 4) != ".BAK"
+      && fieldFile.substr(len - 4) != ".old"))
+      && (len < 5 || fieldFile.substr(len - 5) != ".save"))
       {
-      tempFace.faceIndex = this->FacesOwnerCell->value[i][j];
-      tempFace.neighborFace = false;
-      this->FacesOfCell->value[i].push_back(tempFace);
+      vtkFoamIOobject io;
+      if(io.open(tempPath + "/" + fieldFile)) // file exists and readable
+        {
+        if(io.className() == "volScalarField" ||
+          io.className() == "volVectorField")
+          {
+          // real file name
+          this->TimeStepData->value.push_back(fieldFile);
+          // object name
+          this->CellDataArraySelection->AddArray(io.objectName().c_str());
+          }
+        io.close();
+        }
+#if 0
+      // warning is turned off in favor of silence
+      else
+        {
+        vtkWarningMacro(<< "File " << io.fileName().c_str()
+          << " is not a valid OpenFOAM object. Reason: "
+          << io.error().str().c_str() << " in line "
+          << io.lineNumber());
+        }
+#endif
       }
     }
+  directory->Delete();
 
-  //add neighbor faces to cell
-  for(int i = 0; i < (int)this->FacesNeighborCell->value.size(); i++)
-    {
-    for(int j = 0; j < (int)this->FacesNeighborCell->value[i].size(); j++)
-      {
-      tempFace.faceIndex = this->FacesNeighborCell->value[i][j];
-      tempFace.neighborFace = true;
-      this->FacesOfCell->value[i].push_back(tempFace);
+  // locate laglangian fields
+
+  //clear prior timestep data
+  this->LagrangianTimeStepData->value.clear();
+  // Commented out the "RemoveAllArrays()" line....
+  //this->PointDataArraySelection->RemoveAllArrays();
+
+  //open the directory and get num of files
+  tempPath += "/lagrangian";
+  directory = vtkDirectory::New();
+  if(directory->Open(tempPath.c_str()))
+    {
+    nFieldFiles = directory->GetNumberOfFiles();
+    for(int j = 0; j < nFieldFiles; j++)
+      {
+      const vtkStdString fieldFile(directory->GetFile(j));
+      const int len = fieldFile.length();
+
+      // excluded extensions based on src/OpenFOAM/OSspecific/Unix/Unix.C
+      if(!directory->FileIsDirectory(fieldFile.c_str())
+        && fieldFile.substr(len - 1) != "~"
+        && (len < 4 || (fieldFile.substr(len - 4) != ".bak"
+        && fieldFile.substr(len - 4) != ".BAK"
+        && fieldFile.substr(len - 4) != ".old"))
+        && (len < 5 || fieldFile.substr(len - 5) != ".save"))
+        {
+        vtkFoamIOobject io;
+        if(io.open(tempPath + "/" + fieldFile))
+          {
+          if(io.className() == "scalarField"
+            || io.className() == "vectorField")
+            {
+            // real file name
+            this->LagrangianTimeStepData->value.push_back(fieldFile);
+            // object name
+            this->PointDataArraySelection->AddArray(
+	      (vtkStdString("lagrangian/") + io.objectName()).c_str());
+            }
+          else if(io.className() == "Cloud" && io.objectName() == "positions")
+            {
+            this->PatchDataArraySelection->AddArray("Lagrangian Particles");
+            }
+          io.close();
+          }
+        }
       }
     }
+  directory->Delete();
 
-  //clean up memory
-  this->FacesOwnerCell->value.clear();
-  this->FacesNeighborCell->value.clear();
-  return;
+  // refresh selection status
+  this->SelectionModified();
 }
 
 // ****************************************************************************
-//  Method: vtkOpenFOAMReader::MakeInternalMesh
+//  Method: vtkOpenFOAMReader::InsertCellToGrid
 //
 //  Purpose:
-//  derive cell types and create the internal mesh
+//  determine cell shape and insert the cell into the mesh
+//  hexahedron, prism, pyramid, tetrahedron, wedge&tetWedge
 //
 // ****************************************************************************
-vtkUnstructuredGrid * vtkOpenFOAMReader::MakeInternalMesh()
+void vtkOpenFOAMReader::InsertCellToGrid(vtkUnstructuredGrid* internalMesh,
+  int cellId, faceVectorVector *facesOfCell, intVectorVector* facesPoints)
 {
-  vtkDebugMacro(<<"Make internal mesh");
-  //initialize variables
-  bool foundDup = false;
-  vtkstd::vector< int > cellPoints;
-  vtkstd::vector< int > tempFaces[2];
-  vtkstd::vector< int > firstFace;
-  int pivotPoint = 0;
-  int i, j, k, l, pCount;
-  int faceCount = 0;
-
-  //Create Mesh
-  vtkUnstructuredGrid *  internalMesh = vtkUnstructuredGrid::New();
-  //loop through each cell, derive type and insert it into the mesh
-  //hexahedron, prism, pyramid, tetrahedron, wedge&tetWedge
-  for(i = 0; i < (int)this->FacesOfCell->value.size(); i++)  //each cell
-    {
+  // aliases
+  const vtkstd::vector<vtkstd::vector<int> >& facePoints
+    = facesPoints->value;
+  const vtkstd::vector<face>& cellFaces = facesOfCell->value[cellId];
+  const size_t nCellFaces = cellFaces.size();
 
-    //calculate the total points for the cell
-    //used to derive cell type
-    int totalPointCount = 0;
-    for(j = 0; j < (int)this->FacesOfCell->value[i].size(); j++)  //each face
-      {
-      totalPointCount += 
-        (int)this->FacePoints->
-                 value[this->FacesOfCell->value[i][j].faceIndex].size();
-      }
+  vtkIdList* cellPoints = vtkIdList::New();
 
-    // using cell type - order points, create cell, & add to mesh
-    //OFhex | vtkHexahedron
-    if (totalPointCount == 24)
+  // determine type of the cell
+  // based on src/OpenFOAM/meshes/meshShapes/cellMatcher/{hex|prism|pyr|tet}-
+  // Matcher.C
+  int cellType = VTK_CONVEX_POINT_SET;
+  if(nCellFaces == 6)
+    {
+    size_t j = 0;
+    for(j = 0; j < nCellFaces; j++)
       {
-      faceCount = 0;
-
-      //get first face
-      for(j = 0; j <
-        (int)this->FacePoints->
-                    value[this->FacesOfCell->value[i][0].faceIndex].size(); j++)
+      if(facePoints[cellFaces[j].faceIndex].size() != 4)
         {
-        firstFace.push_back(this->FacePoints->value[
-          this->FacesOfCell->value[i][0].faceIndex][j]);
+        break;
         }
-
-      //patch: if it is a neighbor face flip the points
-      if(this->FacesOfCell->value[i][0].neighborFace)
+      }
+    if(j == nCellFaces)
+      {
+      cellType = VTK_HEXAHEDRON;
+      }
+    }
+  else if(nCellFaces == 5)
+    {
+    int nTris = 0, nQuads = 0;
+    for(size_t j = 0; j < nCellFaces; j++)
+      {
+      const int nPoints = facePoints[cellFaces[j].faceIndex].size();
+      if(nPoints == 3)
         {
-        int tempPop;
-        for(k = 0; k < (int)firstFace.size() - 1; k++)
-          {
-          tempPop = firstFace[firstFace.size()-1];
-          firstFace.pop_back();
-          firstFace.insert(firstFace.begin()+1+k, tempPop);
-          }
+        nTris++;
         }
-
-      //add first face to cell points
-      for(j =0; j < (int)firstFace.size(); j++)
+      else if(nPoints == 4)
         {
-        cellPoints.push_back(firstFace[j]);
+        nQuads++;
         }
-
-      //find the opposite face and order the points correctly
-      for(int pointCount = 0; pointCount < (int)firstFace.size(); pointCount++)
+      else
         {
-
-        //find the other 2 faces containing each point
-        for(j = 1; j < (int)this->FacesOfCell->value[i].size(); j++) //each face
-          {
-          for(k = 0; k < (int)this->FacePoints->value[
-            this->FacesOfCell->value[i][j].faceIndex].size(); k++) //each point
-            {
-            if(firstFace[pointCount] == this->FacePoints->value[
-              this->FacesOfCell->value[i][j].faceIndex][k])
-              {
-              //ANOTHER FACE WITH THE POINT
-              for(l = 0; l < (int)this->FacePoints->value[
-                this->FacesOfCell->value[i][j].faceIndex].size(); l++)
-                {
-                tempFaces[faceCount].push_back(this->FacePoints->value[
-                  this->FacesOfCell->value[i][j].faceIndex][l]);
-                }
-              faceCount++;
-              }
-            }
-          }
-
-        //locate the pivot point contained in faces 0 & 1
-        for(j = 0; j < (int)tempFaces[0].size(); j++)
-          {
-          for(k = 0; k < (int)tempFaces[1].size(); k++)
-            {
-            if(tempFaces[0][j] == tempFaces[1][k] && tempFaces[0][j] !=
-              firstFace[pointCount])
-              {
-              pivotPoint = tempFaces[0][j];
-              break;
-              }
-            }
-          }
-        cellPoints.push_back(pivotPoint);
-        tempFaces[0].clear();
-        tempFaces[1].clear();
-        faceCount=0;
+        break;
         }
-
-      //create the hex cell and insert it into the mesh
-      vtkHexahedron * hexahedron= vtkHexahedron::New();
-      for(pCount = 0; pCount < (int)cellPoints.size(); pCount++)
+      }
+    if(nTris == 2 && nQuads == 3)
+      {
+      cellType = VTK_WEDGE;
+      }
+    else if(nTris == 4 && nQuads == 1)
+      {
+      cellType = VTK_PYRAMID;
+      }
+    }
+  else if(nCellFaces == 4)
+    {
+    size_t j = 0;
+    for(j = 0; j < nCellFaces; j++)
+      {
+      if(facePoints[cellFaces[j].faceIndex].size() != 3)
         {
-        hexahedron->GetPointIds()->SetId(pCount, cellPoints[pCount]);
+        break;
         }
-      internalMesh->InsertNextCell(hexahedron->GetCellType(),
-        hexahedron->GetPointIds());
-      hexahedron->Delete();
-      cellPoints.clear();
-      firstFace.clear();
       }
+    if(j == nCellFaces)
+      {
+      cellType = VTK_TETRA;
+      }
+    }
 
-    //OFprism | vtkWedge
-    else if (totalPointCount == 18)
+  // not an Hex/Wedge/Pyramid/Tetra
+  if(cellType == VTK_CONVEX_POINT_SET)
+    {
+    int nPoints = 0;
+    for(size_t j = 0; j < nCellFaces; j++)
       {
-      faceCount = 0;
-      int index = 0;
+      nPoints += facePoints[cellFaces[j].faceIndex].size();
+      }
+    if(nPoints == 0)
+      {
+      cellType = VTK_EMPTY_CELL;
+      }
+    }
 
-      //find first triangular face
-      for(j = 0; j < (int)this->FacesOfCell->value[i].size(); j++)  //each face
+  // Cell shape constructor based on the one implementd by Terry
+  // Jordan, with some improvements. not as elegant as the one in
+  // OpenFOAM but works reasonably fast.
+
+  // OFhex | vtkHexahedron || OFprism | vtkWedge
+  if (cellType == VTK_HEXAHEDRON || cellType == VTK_WEDGE)
+    {
+    // find the base face number
+    size_t face0Id = 0;
+    if(cellType == VTK_HEXAHEDRON)
+      {
+      cellPoints->SetNumberOfIds(8);
+      face0Id = 0;
+      }
+    else // VTK_WEDGE
+      {
+      cellPoints->SetNumberOfIds(6);
+      for(size_t j = 0; j < nCellFaces; j++)
         {
-        if((int)this->FacePoints->
-            value[this->FacesOfCell->value[i][j].faceIndex].size() == 3)
+        if(facePoints[cellFaces[j].faceIndex].size() == 3)
           {
-          for(k = 0; k < (int)this->FacePoints->value[
-              this->FacesOfCell->value[i][j].faceIndex].size(); k++)
-            {
-            firstFace.push_back(this->FacePoints->value[
-              this->FacesOfCell->value[i][j].faceIndex][k]);
-            index = j;
-            }
+          face0Id = j;
           break;
           }
         }
+      }
 
+    //get first face in correct order
+    const vtkstd::vector<int>& face0Points
+      = facePoints[cellFaces[face0Id].faceIndex];
+    size_t nFace0Points = face0Points.size();
+    vtkstd::vector<int> baseFace(nFace0Points); // ordered point list
+    if(cellFaces[face0Id].neighborFace)
+      {
       //patch: if it is a neighbor face flip the points
-      if(this->FacesOfCell->value[i][0].neighborFace)
+      for(size_t j = 0; j < nFace0Points; j++)
         {
-        int tempPop;
-        for(k = 0; k < (int)firstFace.size() - 1; k++)
-          {
-          tempPop = firstFace[firstFace.size()-1];
-          firstFace.pop_back();
-          firstFace.insert(firstFace.begin()+1+k, tempPop);
-          }
+        baseFace[j] = face0Points[nFace0Points - 1 - j];
+        //add base face to cell points
+        cellPoints->SetId(j, baseFace[j]);
         }
-
-      //add first face to cell points
-      for(j =0; j < (int)firstFace.size(); j++)
+      }
+    else
+      {
+      for(size_t j = 0; j < nFace0Points; j++)
         {
-        cellPoints.push_back(firstFace[j]);
+        baseFace[j] = face0Points[j];
+        //add base face to cell points
+        cellPoints->SetId(j, baseFace[j]);
         }
+      }
 
-      //find the opposite face and order the points correctly
-      for(int pointCount = 0; pointCount < (int)firstFace.size(); pointCount++)
+    //find the opposite face points and order the points correctly
+    const size_t nBaseFacePoints = baseFace.size();
+    for(size_t pointCount = 0; pointCount < nBaseFacePoints; pointCount++)
+      {
+      int pivotPoint = -1;
+      const int baseFacePointI = baseFace[pointCount];
+      for(size_t j = 0; j < nCellFaces; j++) //each face
         {
-        //find the 2 other faces containing each point
-        for(j = 0; j < (int)this->FacesOfCell->value[i].size(); j++) //each face
+        if(j == face0Id)
+          {
+          continue;
+          }
+        const vtkstd::vector<int>& faceJPoints
+          = facePoints[cellFaces[j].faceIndex];
+        const size_t sizeK = faceJPoints.size();
+        if(sizeK == 3) // omit the other triangular wedge face
+          {
+          continue;
+          }
+        for(size_t k = 0; k < sizeK; k++) //each point
           {
-          for(k = 0; k < (int)this->FacePoints->value[
-              this->FacesOfCell->value[i][j].faceIndex].size(); k++)
+          // if the point matches the point of the base face...
+          if(baseFacePointI == faceJPoints[k])
             {
-            if(firstFace[pointCount] == this->FacePoints->value[
-               this->FacesOfCell->value[i][j].faceIndex][k] && j != index)
+            const bool isNeighbor = cellFaces[j].neighborFace;
+            const int faceJPrevPoint = faceJPoints[(sizeK + k - 1) % sizeK];
+            const int faceJNextPoint = faceJPoints[(k + 1) % sizeK];
+            const int baseFacePrevPoint = baseFace[(nBaseFacePoints
+              + pointCount - 1) % nBaseFacePoints];
+            const int baseFaceNextPoint
+              = baseFace[(pointCount + 1) % nBaseFacePoints];
+
+            // if the next point of the j-th face matches the
+            // previous point of the base face use the previous
+            // point of the j-th face as the pivot point and use the
+            // next point otherwise
+            if(faceJNextPoint == (isNeighbor ? baseFaceNextPoint
+              : baseFacePrevPoint))
               {
-              //ANOTHER FACE WITH POINT
-              for(l = 0; l < (int)this->FacePoints->value[
-                  this->FacesOfCell->value[i][j].faceIndex].size(); l++)
-                {
-                tempFaces[faceCount].push_back(this->FacePoints->value[
-                  this->FacesOfCell->value[i][j].faceIndex][l]);
-                }
-              faceCount++;
+              pivotPoint = faceJPrevPoint;
               }
-            }
-          }
+            else
+              {
+              pivotPoint = faceJNextPoint;
+              }
+            cellPoints->SetId(nBaseFacePoints + pointCount, pivotPoint);
 
-        //locate the pivot point of faces 0 & 1
-        for(j = 0; j < (int)tempFaces[0].size(); j++)
-          {
-          for(k = 0; k < (int)tempFaces[1].size(); k++)
-            {
-            if(tempFaces[0][j] == tempFaces[1][k] && tempFaces[0][j] !=
-              firstFace[pointCount])
+            // we further look at one more vertex:
+            // if the j-th face also shares next point of the base face
+            // add the corresponding opposite point to cell points
+            if(pointCount < nBaseFacePoints - 1 && baseFaceNextPoint
+              == (isNeighbor ? faceJNextPoint : faceJPrevPoint))
               {
-              pivotPoint = tempFaces[0][j];
-              break;
+              pointCount++;
+              cellPoints->SetId(nBaseFacePoints + pointCount,
+                faceJPoints[(k + 2) % sizeK]);
               }
+            break; // look no more points as to the face
             }
           }
-        cellPoints.push_back(pivotPoint);
-        tempFaces[0].clear();
-        tempFaces[1].clear();
-        faceCount=0;
-        }
-
-      //create the wedge cell and insert it into the mesh
-      vtkWedge * wedge= vtkWedge::New();
-      for(pCount = 0; pCount < (int)cellPoints.size(); pCount++)
-        {
-        wedge->GetPointIds()->SetId(pCount, cellPoints[pCount]);
-        }
-      internalMesh->InsertNextCell(wedge->GetCellType(),
-        wedge->GetPointIds());
-      cellPoints.clear();
-      wedge->Delete();
-      firstFace.clear();
-      }
-
-    //OFpyramid | vtkPyramid
-    else if (totalPointCount == 16)
-      {
-      foundDup = false;
-
-      //find the quadratic face
-      for(j = 0; j < (int)this->FacesOfCell->value[i].size(); j++)  //each face
-        {
-        if((int)this->FacePoints->
-            value[this->FacesOfCell->value[i][j].faceIndex].size() == 4)
+        if(pivotPoint >= 0) // break if the pivot point is found
           {
-          for(k = 0; k < (int)this->FacePoints->value[
-              this->FacesOfCell->value[i][j].faceIndex].size(); k++)
-            {
-            cellPoints.push_back(this->FacePoints->value[
-              this->FacesOfCell->value[i][j].faceIndex][k]);
-            }
           break;
           }
         }
+      }
+    //create the hex cell and insert it into the mesh
+    internalMesh->InsertNextCell(cellType, cellPoints);
+    }
 
-      //compare first face points to other faces
-      for(j = 0; j < (int)cellPoints.size(); j++) //each point
+  //OFpyramid | vtkPyramid || OFtet | vtkTetrahedron
+  else if (cellType == VTK_PYRAMID || cellType == VTK_TETRA)
+    {
+    int baseFaceId = -1;
+    if(cellType == VTK_PYRAMID)
+      {
+      for(size_t j = 0; j < nCellFaces; j++)
         {
-        for(k = 0; k < (int)this->FacePoints->value[
-          this->FacesOfCell->value[i][1].faceIndex].size(); k++)
-          {
-          if(cellPoints[j] == this->FacePoints->value[
-            this->FacesOfCell->value[i][1].faceIndex][k])
-            {
-            foundDup = true;
-            }
-          }
-        if(!foundDup)
+        if(facePoints[cellFaces[j].faceIndex].size() == 4)
           {
-          cellPoints.push_back(this->FacePoints->value[
-            this->FacesOfCell->value[i][j].faceIndex][k]);
+          baseFaceId = j;
           break;
           }
         }
-
-      //create the pyramid cell and insert it into the mesh
-      vtkPyramid * pyramid = vtkPyramid::New();
-      for(pCount = 0; pCount < (int)cellPoints.size(); pCount++)
-        {
-        pyramid->GetPointIds()->SetId(pCount, cellPoints[pCount]);
-        }
-      internalMesh->InsertNextCell(pyramid->GetCellType(),
-        pyramid->GetPointIds());
-      cellPoints.clear();
-      pyramid->Delete();
+      cellPoints->SetNumberOfIds(5);
       }
-
-    //OFtet | vtkTetrahedron
-    else if (totalPointCount == 12)
+    else // VTK_TETRA
       {
-      foundDup = false;
+      baseFaceId = 0;
+      cellPoints->SetNumberOfIds(4);
+      }
 
-      //add first face to cell points
-      for(j = 0; j < (int)this->FacePoints->value[
-        this->FacesOfCell->value[i][0].faceIndex].size(); j++)
+    //add first face to cell points
+    const vtkstd::vector<int>& baseFacePoints
+      = facePoints[cellFaces[baseFaceId].faceIndex];
+    const size_t nBaseFacePoints = baseFacePoints.size();
+    if(cellFaces[baseFaceId].neighborFace)
+      {
+      // if it is a neighbor face flip the points
+      for(size_t j = 0; j < nBaseFacePoints; j++)
         {
-        cellPoints.push_back(this->FacePoints->value[
-          this->FacesOfCell->value[i][0].faceIndex][j]);
+        cellPoints->SetId(j, baseFacePoints[nBaseFacePoints - 1 - j]);
         }
-
-      //compare first face to the points of second face
-      for(j = 0; j < (int)cellPoints.size(); j++) //each point
+      }
+    else
+      {
+      for(size_t j = 0; j < nBaseFacePoints; j++)
         {
-        for(k = 0; k < (int)this->FacePoints->value[
-          this->FacesOfCell->value[i][1].faceIndex].size(); k++)
-          {
-          if(cellPoints[j] == this->FacePoints->value[
-            this->FacesOfCell->value[i][1].faceIndex][k])
-            {
-            foundDup = true;
-            }
-          }
-        if(!foundDup)
-          {
-          cellPoints.push_back(this->FacePoints->value[
-            this->FacesOfCell->value[i][j].faceIndex][k]);
-          break;
-          }
+        cellPoints->SetId(j, baseFacePoints[j]);
         }
+      }
 
-      //create the wedge cell and insert it into the mesh
-      vtkTetra * tetra = vtkTetra::New();
-      for(pCount = 0; pCount < (int)cellPoints.size(); pCount++)
+    // compare an adjacent face (any non base face is ok) point 1 to
+    // base face points
+    const int adjacentFaceId = (baseFaceId == 0) ? 1 : baseFaceId - 1;
+    const vtkstd::vector<int>& adjacentFacePoints
+      = facePoints[cellFaces[adjacentFaceId].faceIndex];
+    const int adjacentFacePoint1 = adjacentFacePoints[1];
+    bool foundDup = false;
+    for(size_t j = 0; j < nBaseFacePoints; j++)
+      {
+      // if point 1 of the adjacent face matchs point j of the base face...
+      if(cellPoints->GetId(j) == adjacentFacePoint1)
         {
-        tetra->GetPointIds()->SetId(pCount, cellPoints[pCount]);
+        // if point 2 of the adjacent face matches the previous point
+        // of the base face use point 0 of the adjacent face as the
+        // pivot point; use point 2 otherwise
+        cellPoints->SetId(nBaseFacePoints, (adjacentFacePoints[2]
+          == cellPoints->GetId((cellFaces[adjacentFaceId].neighborFace
+          ? (j + 1) : (nBaseFacePoints + j - 1)) % nBaseFacePoints)) ? 
+          adjacentFacePoints[0] : adjacentFacePoints[2]);
+        foundDup = true;
+        break;
         }
-      internalMesh->InsertNextCell(tetra->GetCellType(),
-        tetra->GetPointIds());
-      cellPoints.clear();
-      tetra->Delete();
       }
-
-    //erronous cells
-    else if(totalPointCount == 0)
+    // if point 1 of the adjacent face does not match any points of
+    // the base face, it's the pivot point
+    if(!foundDup)
       {
-      vtkWarningMacro("Warning: No points in cell.");
+      cellPoints->SetId(nBaseFacePoints, adjacentFacePoint1);
       }
 
-    //OFpolyhedron || vtkConvexPointSet
-    else
-      {
-      vtkWarningMacro("Warning: Polyhedral Data is very Slow!");
-      foundDup = false;
+    //create the tetra cell and insert it into the mesh
+    internalMesh->InsertNextCell(cellType, cellPoints);
+    }
 
-      //get first face
-      for(j = 0; j < (int)this->FacePoints->value[
-        this->FacesOfCell->value[i][0].faceIndex].size(); j++)
+  //erronous cells
+  else if(cellType == VTK_EMPTY_CELL)
+    {
+    vtkWarningMacro("Warning: No points in cell.");
+    }
+
+  //OFpolyhedron || vtkConvexPointSet
+  else
+    {
+    // vtkWarningMacro("Warning: Polyhedral Data is very Slow!");
+
+    //get first face
+    const vtkstd::vector<int>& baseFacePoints
+      = facePoints[cellFaces[0].faceIndex];
+    const size_t nBaseFacePoints = baseFacePoints.size();
+    //add first face to cell points
+    // not sure if flipping is necessary but do it anyway
+    if(cellFaces[0].neighborFace)
+      {
+      // if it is a neighbor face flip the points
+      for(size_t j = 0; j < nBaseFacePoints; j++)
         {
-        firstFace.push_back(this->FacePoints->value[
-          this->FacesOfCell->value[i][0].faceIndex][j]);
+        cellPoints->InsertNextId(baseFacePoints[nBaseFacePoints - 1 - j]);
         }
-
-      //add first face to cell points
-      for(j =0; j < (int)firstFace.size(); j++)
+      }
+    else
+      {
+      for(size_t j = 0; j < nBaseFacePoints; j++)
         {
-        cellPoints.push_back(firstFace[j]);
+        cellPoints->InsertNextId(baseFacePoints[j]);
         }
+      }
 
-      //loop through faces and create a list of all points
-      //j = 1 skip firstFace
-      for(j = 1; j < (int)this->FacesOfCell->value[i].size(); j++)
+    //loop through faces and create a list of all points
+    //j = 1 skip baseFace
+    for(size_t j = 1; j < nCellFaces; j++)
+      {
+      //remove duplicate points from faces
+      const vtkstd::vector<int>& faceJPoints
+        = facePoints[cellFaces[j].faceIndex];
+      const size_t nFaceJPoints = faceJPoints.size();
+      for(size_t k = 0; k < nFaceJPoints; k++)
         {
-        //remove duplicate points from faces
-        for(k = 0; k < (int)this->FacePoints->value[
-          this->FacesOfCell->value[i][j].faceIndex].size(); k++)
+        const int faceJPointK = faceJPoints[k];
+        bool foundDup = false;
+        for(int l = 0; l < cellPoints->GetNumberOfIds(); l++)
           {
-          for(l = 0; l < (int)cellPoints.size(); l++);
-            {
-            if(cellPoints[l] == this->FacePoints->value[
-              this->FacesOfCell->value[i][j].faceIndex][k])
-              {
-              foundDup = true;
-              }
-            }
-          if(!foundDup)
+          if(cellPoints->GetId(l) == faceJPointK)
             {
-            cellPoints.push_back(this->FacePoints->value[
-              this->FacesOfCell->value[i][j].faceIndex][k]);
-            foundDup = false;
+            foundDup = true;
+            break; // look no more
             }
           }
+        if(!foundDup)
+          {
+          cellPoints->InsertNextId(faceJPointK);
+          }
         }
-
-      //create the poly cell and insert it into the mesh
-      vtkConvexPointSet * poly = vtkConvexPointSet::New();
-      poly->GetPointIds()->SetNumberOfIds(cellPoints.size());
-      for(pCount = 0; pCount < (int)cellPoints.size(); pCount++)
-        {
-        poly->GetPointIds()->SetId(pCount, cellPoints[pCount]);
-        }
-      internalMesh->InsertNextCell(poly->GetCellType(),
-        poly->GetPointIds());
-      cellPoints.clear();
-      firstFace.clear();
-      poly->Delete();
       }
-    }
 
-  //set the internal mesh points
-  internalMesh->SetPoints(Points);
-  vtkDebugMacro(<<"Internal mesh made");
-  return internalMesh;
+    //create the poly cell and insert it into the mesh
+    internalMesh->InsertNextCell(VTK_CONVEX_POINT_SET, cellPoints);
+    }
+  cellPoints->Delete();
 }
 
 // ****************************************************************************
-//  Method: vtkOpenFOAMReader::ControlDictDataParser
+//  Method: vtkOpenFOAMReader::MakeInternalMesh
 //
 //  Purpose:
-//  parse out double values for controlDict entries
-//  utility function
+//  derive cell types and create the internal mesh
 //
 // ****************************************************************************
-double vtkOpenFOAMReader::ControlDictDataParser(const char * lineIn)
+vtkUnstructuredGrid * vtkOpenFOAMReader::MakeInternalMesh(
+  faceVectorVector *facesOfCell, intVectorVector* facesPoints,
+  vtkPoints* points)
 {
-  double value;
-  vtkstd::string line(lineIn);
-  line.erase(line.begin()+line.find(";"));
-  vtkstd::string token;
-  vtksys_ios::stringstream tokenizer(line);
+  vtkDebugMacro(<<"Make internal mesh");
 
-  //parse to the final entry - double
-  //while(tokenizer>>token);
-  while(!tokenizer.eof())
+  //Create Mesh
+  vtkUnstructuredGrid* internalMesh = vtkUnstructuredGrid::New();
+  internalMesh->Allocate(facesOfCell->value.size());
+
+  //loop through each cell, derive type and insert it into the mesh
+  //hexahedron, prism, pyramid, tetrahedron, wedge&tetWedge
+  for(size_t i = 0; i < facesOfCell->value.size(); i++)  //each cell
     {
-    tokenizer >> token;
+    this->InsertCellToGrid(internalMesh, i, facesOfCell, facesPoints);
     }
+  //set the internal mesh points
+  internalMesh->SetPoints(points);
+  vtkDebugMacro(<<"Internal mesh made");
 
-  vtksys_ios::stringstream conversion(token);
-  conversion >> value;
-  return value;
+  return internalMesh;
 }
 
-// ****************************************************************************
-//  Method: vtkOpenFOAMReader::ReadControlDict
-//
-//  Pupose:
-//  reads the controlDict File
-//  gather the necessary information to create a path to the data
-//
-// ****************************************************************************
-void vtkOpenFOAMReader::ReadControlDict ()
+// List time directories according to controlDict
+bool vtkOpenFOAMReader::ListTimeDirectoriesByControlDict(vtkFoamDict* dictPtr)
 {
-  vtkDebugMacro(<<"Read controlDict");
-  //create variables
-  vtkstd::string temp;
-  double startTime;
-  double endTime;
-  double deltaT;
-  double writeInterval;
-  double timeStepIncrement;
-  vtkstd::string writeControl;
-  vtkstd::string timeFormat;
-  stdString* tempStringStruct;
-
-  ifstream * input = new ifstream(this->Path->value.c_str(), ios::in VTK_IOS_NOCREATE);
-
-  //create the path to the data directory
-  this->PathPrefix->value = this->Path->value;
-  this->PathPrefix->value.erase(this->PathPrefix->value.begin()+
-    this->PathPrefix->value.find("system"),this->PathPrefix->value.end());
-  vtkDebugMacro(<<"Path: "<<this->PathPrefix->value.c_str());
-
-  //find Start Time
-  tempStringStruct = this->GetLine(input);
-  temp = tempStringStruct->value;
-  delete tempStringStruct;
+  vtkFoamDict& dict = *dictPtr;
 
-  //while(!(temp.compare(0,8,"startTime",0,8) == 0))
-  while (strcmp(temp.substr(0,9).c_str(), "startTime"))
+  vtkFoamEntry& startTimeEntry = dict.lookup("startTime");
+  if(!startTimeEntry.found())
     {
-    tempStringStruct = this->GetLine(input);
-    temp = tempStringStruct->value;
-    delete tempStringStruct;
+    vtkErrorMacro(<< "startTime entry not found in controlDict");
+    return false;
     }
-  startTime = this->ControlDictDataParser(temp.c_str());
-  vtkDebugMacro(<<"Start time: "<<startTime);
+  double startTime;
+  startTimeEntry >> startTime;
 
-  //find End Time
-  //while(!(temp.compare(0,6,"endTime",0,6) == 0))
-  while (strcmp(temp.substr(0,7).c_str(), "endTime"))
+  vtkFoamEntry& endTimeEntry = dict.lookup("endTime");
+  if(!endTimeEntry.found())
     {
-    tempStringStruct = this->GetLine(input);
-    temp = tempStringStruct->value;
-    delete tempStringStruct;
+    vtkErrorMacro(<< "endTime entry not found in controlDict");
+    return false;
     }
-  endTime = this->ControlDictDataParser(temp.c_str());
-  vtkDebugMacro(<<"End time: "<<endTime);
+  double endTime;
+  endTimeEntry >> endTime;
 
-  //find Delta T
-  //while(!(temp.compare(0,5,"deltaT",0,5) == 0))
-  while (strcmp(temp.substr(0,6).c_str(), "deltaT"))
+  vtkFoamEntry& deltaTEntry = dict.lookup("deltaT");
+  if(!deltaTEntry.found())
     {
-    tempStringStruct = this->GetLine(input);
-    temp = tempStringStruct->value;
-    delete tempStringStruct;
+    vtkErrorMacro(<< "deltaT entry not found in controlDict");
+    return false;
     }
-  deltaT = this->ControlDictDataParser(temp.c_str());
-  vtkDebugMacro(<<"deltaT: "<<deltaT);
+  double deltaT;
+  deltaTEntry >> deltaT;
 
-  //find write control
-  //while(!(temp.compare(0,11,"writeControl",0,11) == 0))
-  while (strcmp(temp.substr(0,12).c_str(), "writeControl"))
+  vtkFoamEntry& writeIntervalEntry = dict.lookup("writeInterval");
+  if(!writeIntervalEntry.found())
     {
-    tempStringStruct = this->GetLine(input);
-    temp = tempStringStruct->value;
-    delete tempStringStruct;
+    vtkErrorMacro(<< "writeInterval entry not found in controlDict");
+    return false;
     }
+  double writeInterval;
+  writeIntervalEntry >> writeInterval;
 
-  temp.erase(temp.begin()+temp.find(";"));
-  vtkstd::string token;
-  vtksys_ios::stringstream tokenizer(temp);
-
-  //while(tokenizer >> token);
-  while(!tokenizer.eof())
+  vtkFoamEntry& timeFormatEntry = dict.lookup("timeFormat");
+  if(!timeFormatEntry.found())
     {
-    tokenizer >> token;
+    vtkErrorMacro(<< "timeFormat entry not found in controlDict");
+    return false;
     }
-  writeControl = token;
-  vtkDebugMacro(<<"Write control: "<<writeControl.c_str());
+  vtkStdString timeFormat;
+  timeFormatEntry >> timeFormat;
 
-  //find write interval
-  //while(!(temp.compare(0,12,"writeInterval",0,12) == 0))
-  while (strcmp(temp.substr(0,13).c_str(), "writeInterval"))
+  vtkFoamEntry& timePrecisionEntry = dict.lookup("timePrecision");
+  int timePrecision;
+  if(timePrecisionEntry.found())
+    {
+    timePrecisionEntry >> timePrecision;
+    }
+  else
     {
-    tempStringStruct = this->GetLine(input);
-    temp = tempStringStruct->value;
-    delete tempStringStruct;
+    timePrecision = 6; // the default value
     }
-  writeInterval = this->ControlDictDataParser(temp.c_str());
-  vtkDebugMacro(<<"Write interval: "<<writeInterval);
 
   //calculate the time step increment based on type of run
-  //if(writeControl.compare(0,7,"timeStep",0,7) == 0)
-  if(!strcmp(writeControl.substr(0,8).c_str(), "timeStep"))
+  vtkStdString writeControl;
+  dict.lookup("writeControl") >> writeControl;
+  double timeStepIncrement;
+  if(writeControl == "timeStep")
     {
     vtkDebugMacro(<<"Time step type data");
     timeStepIncrement = writeInterval * deltaT;
     }
-  else
+  else if(writeControl == "runTime" || writeControl == "adjustableRunTime")
     {
     vtkDebugMacro(<<"Run time type data");
     timeStepIncrement = writeInterval;
     }
-
-  //find time format
-  while(temp.find("timeFormat") == vtkstd::string::npos)
+  else
     {
-    tempStringStruct = this->GetLine(input);
-    temp = tempStringStruct->value;
-    delete tempStringStruct;
+    vtkErrorMacro(<<"Time step can't be determined because writeControl is"
+      " set to " << writeControl.c_str());
+    return false;
     }
-  timeFormat = temp;
 
   //calculate how many timesteps there should be
-  float tempResult = ((endTime-startTime)/timeStepIncrement);
-  int tempNumTimeSteps = (int)(tempResult+0.5)+1;  //+0.1 to round up
+  double tempResult = (endTime - startTime) / timeStepIncrement;
+  int tempNumTimeSteps = (int)(tempResult + 0.5) + 1;  //+0.1 to round up
+
   //make sure time step dir exists
   vtkstd::vector< double > tempSteps;
   vtkDirectory * test = vtkDirectory::New();
-  vtksys_ios::stringstream parser;
-  double tempStep;
+  this->TimeNames->Initialize();
   for(int i = 0; i < tempNumTimeSteps; i++)
     {
-    tempStep = i*timeStepIncrement + startTime;
-    parser.clear();
-    if(timeFormat.find("general") != vtkstd::string::npos)
+    vtkstd::ostringstream parser;
+    const double tempStep = i * timeStepIncrement + startTime;
+
+    // determine time name based on Foam::Time::timeName()
+    // in src/OpenFOAM/db/Time/Time.C
+#ifdef _MSC_VER
+    bool correctExponent = true;
+#endif
+    if(timeFormat == "general")
+      {
+      // "do not use std:: or vtkstd:: when using anything declared in
+      // iostream," according to VTK coding standards, but following
+      // the instruction causes errors...
+      parser.setf(vtkstd::ios_base::fmtflags(0), vtkstd::ios_base::floatfield);
+      }
+    else if(timeFormat == "fixed")
+      {
+      parser.setf(vtkstd::ios_base::fmtflags(vtkstd::ios_base::fixed),
+        vtkstd::ios_base::floatfield);
+#ifdef _MSC_VER
+      correctExponent = false;
+#endif
+      }
+    else if(timeFormat == "scientific")
+      {
+      parser.setf(vtkstd::ios_base::fmtflags(vtkstd::ios_base::scientific),
+        vtkstd::ios_base::floatfield);
+      }
+    else
       {
-      parser << tempStep;
+      vtkWarningMacro("Warning: unsupported time format. Assuming general.");
+      parser.setf(vtkstd::ios_base::fmtflags(0), vtkstd::ios_base::floatfield);
+      }
+    parser.precision(timePrecision);
+    parser << tempStep; // stringstream doesn't require ends
+#ifdef _MSC_VER
+    // workaround for format difference in MSVC++:
+    // remove an extra 0 from exponent
+    if(correctExponent)
+      {
+      vtkStdString tempStr(parser.str());
+      size_t pos = tempStr.find('e');
+      if(pos != vtkStdString::npos && tempStr.length() >= pos + 3
+        && tempStr[pos + 2] == '0')
+        {
+        tempStr.erase(pos + 2, 1);
+        parser.str(tempStr);
+        }
       }
-    else
+#endif
+    if(test->Open((*this->PathPrefix + parser.str()).c_str()))
       {
-      parser << ios::scientific <<tempStep;
+      tempSteps.push_back(tempStep);
+      this->TimeNames->InsertNextValue(parser.str());
       }
-    if(test->Open((this->PathPrefix->value+parser.str().c_str()).c_str()))
+    // necessary for reading the case/0 directory whatever the timeFormat is
+    // based on Foam::Time::operator++() in src/OpenFOAM/db/Time/Time.C
+    else if((fabs(tempStep) < 1.0e-14L) // 10*SMALL
+      && test->Open((*this->PathPrefix + vtkStdString("0")).c_str()))
       {
       tempSteps.push_back(tempStep);
+      this->TimeNames->InsertNextValue(vtkStdString("0"));
       }
     }
   test->Delete();
+  this->TimeNames->Squeeze();
 
   //Add the time steps that actually exist to steps
   //allows the run to be stopped short of controlDict spec
   //allows for removal of timesteps
   this->NumberOfTimeSteps = tempSteps.size();
-  this->Steps = new double[this->NumberOfTimeSteps];
-  for(int i = 0; i < this->NumberOfTimeSteps; i++)
+  if(this->NumberOfTimeSteps > 0)
+    {
+    this->Steps = new double[this->NumberOfTimeSteps];
+    for(int i = 0; i < this->NumberOfTimeSteps; i++)
+      {
+      this->Steps[i] = tempSteps[i];
+      }
+    }
+  else
     {
-    this->Steps[i] =tempSteps[i];
+    // dummy for safely deleting later
+    this->Steps = new double[1];
+    this->Steps[0] = 0.0;
     }
+  return true;
+}
 
-  input->close();
-  delete input;
-  vtkDebugMacro(<<"controlDict read");
-  return;
+// list time directories by searching all valid time instances in a
+// case directory
+bool vtkOpenFOAMReader::ListTimeDirectoriesByInstances()
+{
+  // open the case directory
+  vtkDirectory* test = vtkDirectory::New();
+  if(!test->Open(this->PathPrefix->c_str()))
+    {
+    test->Delete();
+    vtkErrorMacro(<<"Can't open directory "<<this->PathPrefix->c_str());
+    return false;
+    }
+
+  // search all the directories in the case directory and detect
+  // directories with names convertible to numbers
+  this->TimeNames->Initialize();
+  vtkDoubleArray* timeValues = vtkDoubleArray::New();
+  const int nFiles = test->GetNumberOfFiles();
+  for(int i = 0; i < nFiles; i++)
+    {
+    const vtkStdString dir = test->GetFile(i);
+    if(test->FileIsDirectory(dir.c_str()))
+      {
+      // check if the name is convertible to a number
+      bool isTimeDir = true;
+      for(size_t j = 0; j < dir.length(); j++)
+        {
+        const char c = dir[j];
+        if(!isdigit(c) && c != '+' && c != '-' && c != '.' && c != 'e'
+          && c != 'E')
+          {
+          isTimeDir = false;
+          break;
+          }
+        }
+      if(!isTimeDir)
+        {
+        continue;
+        }
+
+      // convert to a number
+      char *endptr;
+      double timeValue = strtod(dir.c_str(), &endptr);
+      // check if the value really was converted to a number
+      if(timeValue == 0.0 && endptr == dir.c_str())
+        {
+        continue;
+        }
+
+      // add to the instance list
+      timeValues->InsertNextValue(timeValue);
+      this->TimeNames->InsertNextValue(dir);
+      }
+    }
+  test->Delete();
+  this->TimeNames->Squeeze();
+
+  // sort the detected time directories and set as timesteps
+  this->NumberOfTimeSteps = this->TimeNames->GetSize();
+  if(this->NumberOfTimeSteps > 0)
+    {
+    if(this->NumberOfTimeSteps > 1)
+      {
+      vtkSortDataArray::Sort(timeValues, this->TimeNames);
+      }
+    this->Steps = new double[this->NumberOfTimeSteps];
+    for(int i = 0; i < this->NumberOfTimeSteps; i++)
+      {
+      this->Steps[i] = timeValues->GetValue(i);
+      }
+    }
+  else
+    {
+    // dummy for safely deleting later
+    this->Steps = new double[1];
+    this->Steps[0] = 0.0;
+    }
+  timeValues->Delete();
+
+  return true;
 }
 
 // ****************************************************************************
-//  Method: vtkOpenFOAMReader::GetPoints
+//  Method: vtkOpenFOAMReader::ReadControlDict
 //
-//  Purpose:
-//  read the points file into a vtkPoints
+//  Pupose:
+//  reads the controlDict File
+//  gather the necessary information to create a path to the data
 //
 // ****************************************************************************
-void vtkOpenFOAMReader::GetPoints (int timeState)
+bool vtkOpenFOAMReader::ReadControlDict(const char* pathIn)
 {
-  //path to points file
-  vtkstd::string pointPath = this->PathPrefix->value +
-                             this->PolyMeshPointsDir->value[timeState] +
-                             "/polyMesh/points";
-  vtkDebugMacro(<<"Read points file: "<<pointPath.c_str());
+  vtkStdString path(pathIn);
+  vtkFoamIOobject io;
 
-  vtkstd::string temp;
-  bool binaryWriteFormat;
-  stdString* tempStringStruct;
-  ifstream * input = new ifstream(pointPath.c_str(), ios::in VTK_IOS_NOCREATE);
-  //make sure file exists
-  if(input->fail())
+  if(!io.open(path))
     {
-    input->close();
-    delete input;
-    return;
+    vtkErrorMacro(<<"Error opening " << io.fileName().c_str() << ": "
+      << io.error().str().c_str());
+    return false;
     }
-
-  //determine if file is binary or ascii
-  while(temp.find("format") == vtkstd::string::npos)
+  vtkFoamDict dict;
+  if(!dict.read(io))
+    {
+    vtkErrorMacro(<<"Error reading line " << io.lineNumber()
+      << " of " << io.fileName().c_str() << ": " << io.error().str().c_str());
+    return false;
+    }
+  if(dict.type() != vtkFoamDict::DICTIONARY)
     {
-    tempStringStruct = this->GetLine(input);
-    temp = tempStringStruct->value;
-    delete tempStringStruct;
+    vtkErrorMacro(<<"The file type of " << io.fileName().c_str()
+      << " is not a dictionary");
+    return false;
     }
-  input->close();
 
-  //reopen file in correct format
-  if(temp.find("binary") != vtkstd::string::npos)
+  //create the path to the data directory
+  // note that if a file is specified from command line (with
+  // --data=opt) its path may be a relative one
+  vtkStdString prefix(path);
+  // remove trailing "/controlDict.foam"
+#if defined(_WIN32)
+  size_t pos = prefix.find_last_of("/\\");
+#else
+  size_t pos = prefix.find_last_of('/');
+#endif
+  if(pos != vtkStdString::npos)
     {
-#ifdef _WIN32
-    input->open(pointPath.c_str(), ios::binary | ios::in VTK_IOS_NOCREATE);
+    prefix.erase(pos);
+    // remove trailing "/system" (or any other directory name)
+#if defined(_WIN32)
+    pos = prefix.find_last_of("/\\");
+#else
+    pos = prefix.find_last_of('/');
+#endif
+    if(pos != vtkStdString::npos)
+      {
+      prefix.erase(pos + 1); // preserve the last "/"
+      }
+    else
+      {
+#if defined(_WIN32) && PARAVIEW_VERSION_MAJOR >= 3
+      prefix = ".\\";
 #else
-    input->open(pointPath.c_str(), ios::in VTK_IOS_NOCREATE);
+      prefix = "./";
 #endif
-    binaryWriteFormat = true;
+      }
     }
   else
     {
-    input->open(pointPath.c_str(),ios::in);
-    binaryWriteFormat = false;
-    }
-
-  double x,y,z;
-  vtksys_ios::stringstream tokenizer;
-
-  //instantiate the points class
-  this->Points->Reset();
-
-  //find end of header
-  //while(temp.compare(0,4, "// *", 0, 4) != 0)
-  while (strcmp(temp.substr(0,4).c_str(), "// *"))
-    {
-    tempStringStruct = this->GetLine(input);
-    temp = tempStringStruct->value;
-    delete tempStringStruct;
+#if defined(_WIN32) && PARAVIEW_VERSION_MAJOR >= 3
+    prefix = "..\\";
+#else
+    prefix = "../";
+#endif
     }
+  *this->PathPrefix = prefix;
+  vtkDebugMacro(<<"Path: "<<this->PathPrefix->c_str());
 
-  //find number of points
-  tempStringStruct = this->GetLine(input);
-  temp = tempStringStruct->value;
-  delete tempStringStruct;
-  while(temp.empty())
+  vtkFoamEntry& writeControlEntry = dict.lookup("writeControl");
+  if(!writeControlEntry.found())
     {
-    tempStringStruct = this->GetLine(input);
-    temp = tempStringStruct->value;
-    delete tempStringStruct;
+    vtkErrorMacro(<< "writeControl entry not found in " << pathIn);
+    return false;
     }
-
-  //read number of points
-  tokenizer << temp;
-  tokenizer >> NumPoints;
-  //binary data
-  if(binaryWriteFormat)
+  vtkStdString writeControl;
+  writeControlEntry >> writeControl;
+
+  vtkStdString adjustTimeStep;
+  dict.lookup("adjustTimeStep") >> adjustTimeStep; // empty if not found
+
+  // list time directories according to controlDict if (adjustTimeStep
+  // writeControl) == (off, timeStep) or (on, adjustableRunTime); list
+  // by time instances in the case directory otherwise (different behavior
+  // from paraFoam)
+  // valid switching words taken from src/OpenFOAM/db/Switch/Switch.C
+  if(((adjustTimeStep == "off" || adjustTimeStep == "no"
+    || adjustTimeStep == "n" || adjustTimeStep == "false"
+    || adjustTimeStep == "") && writeControl == "timeStep")
+    || ((adjustTimeStep == "on" || adjustTimeStep == "yes"
+    || adjustTimeStep == "y" || adjustTimeStep == "true")
+    && writeControl == "adjustableRunTime"))
     {
-    input->get(); //parenthesis
-    for(int i = 0; i < NumPoints; i++)
-      {
-      input->read((char *)&x,sizeof(double));
-      input->read((char *)&y,sizeof(double));
-      input->read((char *)&z,sizeof(double));
-      this->Points->InsertPoint(i,x,y,z);
-      }
+    return this->ListTimeDirectoriesByControlDict(&dict);
     }
-
-  //ascii data
   else
     {
-    tempStringStruct = this->GetLine(input);
-    temp = tempStringStruct->value;
-    delete tempStringStruct;
-    for(int i = 0; i < this->NumPoints; i++)
-      {
-      tempStringStruct = this->GetLine(input);
-      temp = tempStringStruct->value;
-      delete tempStringStruct;
-      temp.erase(temp.begin()+temp.find("("));
-      temp.erase(temp.begin()+temp.find(")"));
-      tokenizer.clear();
-      tokenizer << temp;
-      tokenizer >> x;
-      tokenizer >> y;
-      tokenizer >> z;
-      this->Points->InsertPoint(i,x,y,z);
-      }
+    return this->ListTimeDirectoriesByInstances();
     }
-
-  input->close();
-  delete input;
-  vtkDebugMacro(<<"Point file read");
-  return;
 }
 
 // ****************************************************************************
-//  Method: vtkOpenFOAMReader::ReadFacesFile
+//  Method: vtkOpenFOAMReader::GetPoints
 //
 //  Purpose:
-//  read the faces into a vtkstd::vector
+//  read the points file into a vtkPoints
 //
 // ****************************************************************************
-void vtkOpenFOAMReader::ReadFacesFile (const char * facePathIn)
+vtkPoints* vtkOpenFOAMReader::GetPoints (int timeState)
 {
-  vtkstd::string facePath(facePathIn);
-  vtkDebugMacro(<<"Read faces file: "<<facePath.c_str());
-  vtkstd::string temp;
-  stdString* tempStringStruct;
-  bool binaryWriteFormat;
-  ifstream * input = new ifstream(facePath.c_str(), ios::in VTK_IOS_NOCREATE);
-  //make sure file exists
-  if(input->fail())
-    {
-    input->close();
-    delete input;
-    return;
-    }
+  //path to points file
+  vtkStdString pointPath = *this->PathPrefix +
+    this->PolyMeshPointsDir->value[timeState] + "/polyMesh/points";
+  vtkDebugMacro(<<"Read points file: "<<pointPath.c_str());
 
-  //determine if file is binary or ascii
-  while(temp.find("format") == vtkstd::string::npos)
+  vtkFoamIOobject io;
+  if(!(io.open(pointPath) || io.open(pointPath + ".gz")))
     {
-    tempStringStruct = this->GetLine(input);
-    temp = tempStringStruct->value;
-    delete tempStringStruct;
+    vtkErrorMacro(<<"Error opening " << io.fileName().c_str() << ": "
+      << io.error().str().c_str());
+    return NULL;
     }
-  input->close();
-
-  //reopen file in correct format
-  if(temp.find("binary") != vtkstd::string::npos)
+  vtkFoamDict dict;
+  if(!dict.read(io))
     {
-#ifdef _WIN32
-    input->open(facePath.c_str(), ios::binary | ios::in VTK_IOS_NOCREATE);
-#else
-    input->open(facePath.c_str(), ios::in VTK_IOS_NOCREATE);
-#endif
-    binaryWriteFormat = true;
+    vtkErrorMacro(<<"Error reading line " << io.lineNumber()
+      << " of " << io.fileName().c_str() << ": " << io.error().str().c_str());
+    return NULL;
     }
-  else
+  if(dict.type() != vtkFoamDict::VECTORLIST)
     {
-    input->open(facePath.c_str(),ios::in);
-    binaryWriteFormat = false;
+    vtkErrorMacro(<<"The file type of " << io.fileName().c_str()
+      << " is not a vectorList");
+    return NULL;
     }
 
-  vtksys_ios::stringstream tokenizer;
-  size_t pos;
-  int numFacePoints;
-  this->FacePoints->value.clear();
+  vtkDoubleArray& vl = dict.vectorList();
 
-  //find end of header
-  //while(temp.compare(0,4, "// *", 0, 4) != 0)
-  while (strcmp(temp.substr(0,4).c_str(), "// *"))
-    {
-    tempStringStruct = this->GetLine(input);
-    temp = tempStringStruct->value;
-    delete tempStringStruct;
-    }
+  //instantiate the points class
+  vtkPoints* points = vtkPoints::New();
+  const int nPoints = vl.GetNumberOfTuples();
+  points->SetNumberOfPoints(nPoints);
 
-  //find number of faces
-  tempStringStruct = this->GetLine(input);
-  temp = tempStringStruct->value;
-  delete tempStringStruct;
-  while(temp.empty())
+  for(int i = 0, index = 0; i < nPoints; i++, index += 3)
     {
-    tempStringStruct = this->GetLine(input);
-    temp = tempStringStruct->value;
-    delete tempStringStruct;
+    points->SetPoint(i, vl.GetPointer(index));
     }
+  vtkDebugMacro(<<"Point file read");
+  return points;
+}
 
-  //read number of faces
-  tokenizer << temp;
-  tokenizer >> this->NumFaces;
-  this->FacePoints->value.resize(this->NumFaces);
-
-  tempStringStruct = this->GetLine(input);
-  temp = tempStringStruct->value;
-  delete tempStringStruct;//THROW OUT "("
+// ****************************************************************************
+//  Method: vtkOpenFOAMReader::ReadFacesFile
+//
+//  Purpose:
+//  read the faces into a vtkstd::vector
+//
+// ****************************************************************************
+intVectorVector* vtkOpenFOAMReader::ReadFacesFile (const char* facePathIn)
+{
+  const vtkStdString facePath(facePathIn);
+  vtkDebugMacro(<<"Read faces file: "<<facePath.c_str());
 
-  //binary data
-  if(binaryWriteFormat)
+  vtkFoamIOobject io;
+  if(!(io.open(facePath) || io.open(facePath + ".gz")))
     {
-    //char paren;
-    int tempPoint;
-    for(int i = 0; i < NumFaces; i++)
-      {
-      tempStringStruct = this->GetLine(input);
-      temp = tempStringStruct->value;
-      delete tempStringStruct;//THROW OUT blankline
-
-      tempStringStruct = this->GetLine(input);
-      temp = tempStringStruct->value;
-      delete tempStringStruct; //grab point count
-
-      tokenizer.clear();
-      tokenizer << temp;
-      tokenizer >> numFacePoints;
-      this->FacePoints->value[i].resize(numFacePoints);
-      //paren = input->get();  //grab (
-      input->get();  //grab (
-      for(int j = 0; j < numFacePoints; j++)
-        {
-        input->read((char *) &tempPoint, sizeof(int));
-        this->FacePoints->value[i][j] = tempPoint;
-        }
-        tempStringStruct = this->GetLine(input);
-        temp = tempStringStruct->value;
-        delete tempStringStruct; //throw out ) and rest of line
-      }
+    vtkErrorMacro(<<"Error opening " << io.fileName().c_str() << ": "
+      << io.error().str().c_str());
+    return NULL;
     }
-
-  //ascii data
-  else
+  vtkFoamDict dict;
+  if(!dict.read(io))
     {
-    //create vtkstd::vector of points in each face
-    for(int i = 0; i < this->NumFaces; i++)
-      {
-      tempStringStruct = this->GetLine(input);
-      temp = tempStringStruct->value;
-      delete tempStringStruct;
-
-      pos = temp.find("(");
-      vtksys_ios::stringstream ascTokenizer;
-      ascTokenizer << temp.substr(0, pos);
-      temp.erase(0, pos+1);
-      ascTokenizer >> numFacePoints;
-      this->FacePoints->value[i].resize(numFacePoints);
-      for(int j = 0; j < numFacePoints; j++)
-        {
-        pos = temp.find(" ");
-        vtksys_ios::stringstream lineTokenizer;
-        lineTokenizer << temp.substr(0, pos);
-        temp.erase(0, pos+1);
-        lineTokenizer >> this->FacePoints->value[i][j];
-        }
-      }
+    vtkErrorMacro(<<"Error reading line " << io.lineNumber()
+      << " of " << io.fileName().c_str() << ": " << io.error().str().c_str());
+    return NULL;
+    }
+  if(dict.type() != vtkFoamDict::LABELLISTLIST)
+    {
+    vtkErrorMacro(<<"The file type of " << io.fileName().c_str()
+      << " is not a labelListList");
+    return NULL;
     }
 
-  input->close();
-  delete input;
   vtkDebugMacro(<<"Faces read");
-  return;
+  return (intVectorVector *)dict.ptr();
 }
 
 // ****************************************************************************
-//  Method: vtkOpenFOAMReader::ReadOwnerFile
+//  Method: vtkOpenFOAMReader::ReadOwnerNeighborFiles
 //
 //  Purpose:
 //  read the owner file into a vtkstd::vector
 //
 // ****************************************************************************
-void vtkOpenFOAMReader::ReadOwnerFile(const char * ownerPathIn)
+bool vtkOpenFOAMReader::ReadOwnerNeighborFiles(faceVectorVector* facesOfCell,
+  const char* ownerPathIn, const char* neighborPathIn)
 {
-  vtkstd::string ownerPath(ownerPathIn);
+  const vtkStdString ownerPath(ownerPathIn);
   vtkDebugMacro(<<"Read owner file: "<<ownerPath.c_str());
-  vtkstd::string temp;
-  stdString* tempStringStruct;
-  bool binaryWriteFormat;
-  ifstream * input = new ifstream(ownerPath.c_str(), ios::in VTK_IOS_NOCREATE);
-  //make sure file exists
-  if(input->fail())
-    {
-    input->close();
-    delete input;
-    return;
-    }
 
-  //determine if file is binary or ascii
-  while(temp.find("format") == vtkstd::string::npos)
+  vtkFoamIOobject io;
+  if(!(io.open(ownerPath) || io.open(ownerPath + ".gz")))
     {
-    tempStringStruct = this->GetLine(input);
-    temp = tempStringStruct->value;
-    delete tempStringStruct;
+    vtkErrorMacro(<<"Error opening " << io.fileName().c_str() << ": "
+      << io.error().str().c_str());
+    return false;
     }
-  input->close();
-
-  //reopen file in correct format
-  if(temp.find("binary") != vtkstd::string::npos)
+  vtkFoamDict ownerDict;
+  if(!ownerDict.read(io))
     {
-#ifdef _WIN32
-    input->open(ownerPath.c_str(), ios::binary | ios::in VTK_IOS_NOCREATE);
-#else
-    input->open(ownerPath.c_str(), ios::in VTK_IOS_NOCREATE);
-#endif
-    binaryWriteFormat = true;
+    vtkErrorMacro(<<"Error reading line " << io.lineNumber()
+      << " of " << io.fileName().c_str() << ": " << io.error().str().c_str());
+    return false;
     }
-  else
+  if(ownerDict.type() != vtkFoamDict::LABELLIST)
     {
-    input->open(ownerPath.c_str(),ios::in);
-    binaryWriteFormat = false;
+    vtkErrorMacro(<<"The file type of " << io.fileName().c_str()
+      << " is not a labelList");
+    return false;
     }
+  io.close();
 
-  vtkstd::string numFacesStr;
-  int faceValue;
-
-  this->FaceOwner = vtkIntArray::New();
+  const vtkStdString neighborPath(neighborPathIn);
+  vtkDebugMacro(<<"Read neighbor file: "<<neighborPath.c_str());
 
-  vtksys_ios::stringstream tokenizer;
-  tokenizer << this->NumFaces;
-  tokenizer >> numFacesStr;
-  //find end of header & number of faces
-  //while(temp.compare(0,numFacesStr.size(),numFacesStr,
-  //0,numFacesStr.size())!=0)
-  while(strcmp(temp.substr(0,numFacesStr.size()).c_str(), numFacesStr.c_str()))
+  if(!(io.open(neighborPath) || io.open(neighborPath + ".gz")))
     {
-    tempStringStruct = this->GetLine(input);
-    temp = tempStringStruct->value;
-    delete tempStringStruct;
+    vtkErrorMacro(<<"Error opening " << io.fileName().c_str() << ": "
+      << io.error().str().c_str());
+    return false;
     }
-
-  this->FaceOwner->SetNumberOfValues(this->NumFaces);
-
-  //binary data
-  if(binaryWriteFormat)
+  vtkFoamDict neighborDict;
+  if(!neighborDict.read(io))
     {
-    input->get(); //parenthesis
-    for(int i = 0; i < NumFaces; i++)
-      {
-      input->read((char *) &faceValue, sizeof(int));
-      this->FaceOwner->SetValue(i, faceValue);
-      }
+    vtkErrorMacro(<<"Error reading line " << io.lineNumber()
+      << " of " << io.fileName().c_str() << ": " << io.error().str().c_str());
+    return false;
     }
-
-  //ascii data
-  else
+  if(neighborDict.type() != vtkFoamDict::LABELLIST)
     {
-    tempStringStruct = this->GetLine(input);
-    temp = tempStringStruct->value;
-    delete tempStringStruct;
-
-    //read face owners into int array
-    for(int i = 0; i < this->NumFaces; i++)
-      {
-      tempStringStruct = this->GetLine(input);
-      temp = tempStringStruct->value;
-      delete tempStringStruct;
-
-      tokenizer.clear();
-      tokenizer << temp;
-      tokenizer >> faceValue;
-      this->FaceOwner->SetValue(i, faceValue);
-      }
+    vtkErrorMacro(<<"The file type of " << io.fileName().c_str()
+      << " is not a labelList");
+    return false;
     }
 
-  //find the number of cells
-  double  * range;
-  range = this->FaceOwner->GetRange();
-  this->NumCells = (int)range[1]+1;
+  const vtkstd::vector<int>& faceOwner = ownerDict.labelList();
+  const vtkstd::vector<int>& faceNeighbor = neighborDict.labelList();
 
-  //add the face number to the correct cell
-  //according to owner
-  this->FacesOwnerCell->value.resize(this->NumCells);
-  int tempCellId;
-  for(int j = 0; j < this->NumFaces; j++)
+  if(faceOwner.size() != faceNeighbor.size())
     {
-    tempCellId = this->FaceOwner->GetValue(j);
-    if(tempCellId != -1)
-      {
-      this->FacesOwnerCell->value[tempCellId].push_back(j);
-      }
+    vtkErrorMacro(<<"Numbers of faces in owner and neighbor don't match");
+    return false;
     }
 
-  input->close();
-  delete input;
-  vtkDebugMacro(<<"Owner file read");
-  return;
-}
+  const size_t nFaces = faceOwner.size();
 
-// ****************************************************************************
-//  Method: vtkOpenFOAMReader::ReadNeighborFile
-//
-//  Purpose:
-//  read the neighbor file into a vtkstd::vector
-//
-// ****************************************************************************
-void vtkOpenFOAMReader::ReadNeighborFile(const char * neighborPathIn)
-{
-  vtkstd::string neighborPath(neighborPathIn);
-  vtkDebugMacro(<<"Read neighbor file: "<<neighborPath.c_str());
-  vtkstd::string temp;
-  stdString* tempStringStruct;
-  bool binaryWriteFormat;
-  ifstream * input = new ifstream(neighborPath.c_str(), ios::in VTK_IOS_NOCREATE);
-  //make sure file exists
-  if(input->fail())
+  if(nFaces == 0)
     {
-    input->close();
-    delete input;
-    return;
+    vtkWarningMacro(<<"The mesh contains no faces");
     }
 
-  //determine if file is binary or ascii
-  while(temp.find("format") == vtkstd::string::npos)
-    {
-    tempStringStruct = this->GetLine(input);
-    temp = tempStringStruct->value;
-    delete tempStringStruct;
-    }
-  input->close();
+  this->FaceOwner = vtkIntArray::New();
+  vtkIntArray& fo = *this->FaceOwner;
+  fo.SetNumberOfValues(nFaces);
 
-  //reopen file in correct format
-  if(temp.find("binary") != vtkstd::string::npos)
-    {
-#ifdef _WIN32
-    input->open(neighborPath.c_str(), ios::binary | ios::in VTK_IOS_NOCREATE);
-#else
-    input->open(neighborPath.c_str(), ios::in VTK_IOS_NOCREATE);
-#endif
-    binaryWriteFormat = true;
-    }
-  else
+  //read face owners into int array
+  for(size_t i = 0; i < nFaces; i++)
     {
-    input->open(neighborPath.c_str(),ios::in);
-    binaryWriteFormat = false;
-    }
-
-  vtkstd::string numFacesStr;
-  int faceValue;
-  vtkIntArray * faceNeighbor = vtkIntArray::New();
-
-  vtksys_ios::stringstream tokenizer;
-  tokenizer << this->NumFaces;
-  tokenizer >> numFacesStr;
-  //find end of header & number of faces
-  //while(temp.compare(0,numFacesStr.size(),numFacesStr,0,
-  //numFacesStr.size())!=0)
-  while (strcmp(temp.substr(0,numFacesStr.size()).c_str(), numFacesStr.c_str()))
-    {
-    tempStringStruct = this->GetLine(input);
-    temp = tempStringStruct->value;
-    delete tempStringStruct;
+    fo.SetValue(i, faceOwner[i]);
     }
 
-  //read face owners into int array
-  faceNeighbor->SetNumberOfValues(this->NumFaces);
+  //add the face numbers to the correct cell
+  // based on Terry's code and
+  // src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCells.C
 
-  //binary data
-  if(binaryWriteFormat)
+  //find the number of cells
+  int nCells = -1;
+  for(size_t i = 0; i < nFaces; i++)
     {
-    input->get(); //parenthesis
-    for(int i = 0; i < this->NumFaces; i++)
+    if(nCells < faceOwner[i]) // max(nCells, faceOwner[i])
+      {
+      nCells = faceOwner[i];
+      }
+    if(nCells < faceNeighbor[i]) // max(nCells, faceNeighbor[i])
       {
-      input->read((char *) &faceValue, sizeof(int));
-      faceNeighbor->SetValue(i, faceValue);
+      nCells = faceNeighbor[i];
       }
     }
+  nCells++;
 
-  //ascii data
-  else
+  if(nCells == 0)
     {
-    tempStringStruct = this->GetLine(input);
-    temp = tempStringStruct->value;
-    delete tempStringStruct;//throw away (
-    //read face owners into int array
-    for(int i = 0; i < this->NumFaces; i++)
-      {
-      tempStringStruct = this->GetLine(input);
-      temp = tempStringStruct->value;
-      delete tempStringStruct;
-
-      tokenizer.clear();
-      tokenizer << temp;
-      tokenizer >> faceValue;
-      faceNeighbor->SetValue(i, faceValue);
-      }
+    vtkWarningMacro(<<"The mesh contains no cells");
     }
-  //No need to recalulate the Number of Cells
-  this->FacesNeighborCell->value.resize(this->NumCells);
 
-  //add face number to correct cell
-  int tempCellId;
-  for(int j = 0; j < this->NumFaces; j++)
+  this->NumCells = nCells;
+
+  // count number of faces per cell
+  vtkstd::vector<int> nCellFaces(nCells, 0);
+  for(size_t i = 0; i < nFaces; i++)
     {
-    tempCellId = faceNeighbor->GetValue(j);
-    if(tempCellId != -1)
+    // simpleFoam/pitzDaily3Blocks has faces with owner cell number -1
+    if(faceOwner[i] >= 0)
+      {
+      nCellFaces[faceOwner[i]]++;
+      }
+    if(faceNeighbor[i] >= 0)
       {
-      this->FacesNeighborCell->value[tempCellId].push_back(j);
+      nCellFaces[faceNeighbor[i]]++;
       }
     }
 
-  faceNeighbor->Delete();
-  input->close();
-  delete input;
-  vtkDebugMacro(<<"Neighbor file read");
-  return;
-}
-
-// ****************************************************************************
-//  Method: vtkOpenFOAMReader::PopulatePolyMeshDirArrays
-//
-//  Purpose:
-//  create a Lookup Table containing the location of the points
-//  and faces files for each time steps mesh
-//
-// ****************************************************************************
-void vtkOpenFOAMReader::PopulatePolyMeshDirArrays()
-{
-  vtkDebugMacro(<<"Create list of points/faces file directories");
-  vtksys_ios::stringstream path;
-  vtksys_ios::stringstream timeStep;
-  bool facesFound;
-  bool pointsFound;
-  bool polyMeshFound;
-
-  //intialize size to number of timesteps
-  this->PolyMeshPointsDir->value.resize(this->NumberOfTimeSteps);
-  this->PolyMeshFacesDir->value.resize(this->NumberOfTimeSteps);
-
-  //loop through each timestep
-  for(int i = 0; i < this->NumberOfTimeSteps; i++)
+  // size cellFaces
+  vtkstd::vector<vtkstd::vector<face> >& cellFaces = facesOfCell->value;
+  cellFaces.resize(nCells);
+  for(int i = 0; i < nCells; i++)
     {
-    polyMeshFound = false;
-    facesFound = false;
-    pointsFound = false;
-
-    //create the path to the timestep
-    path.clear();
-    timeStep.clear();
-    timeStep << Steps[i];
-    path << this->PathPrefix->value <<timeStep.str() << "/";
+    cellFaces[i].resize(nCellFaces[i]);
+    }
 
-    //get the number of files
-    vtkDirectory * directory = vtkDirectory::New();
-    directory->Open(path.str().c_str());
-    int numFiles = directory->GetNumberOfFiles();
-    //Look for polyMesh Dir
-    for(int j = 0; j < numFiles; j++)
+  //add face numbers to correct cell
+  for(size_t i = 0; i < nFaces; i++)
+    {
+    const int ownerCell = faceOwner[i]; // must be a signed int
+    // simpleFoam/pitzDaily3Blocks has faces with owner cell number -1
+    if(ownerCell >= 0)
       {
-      vtkstd::string tempFile(directory->GetFile(j));
-      if(tempFile.find("polyMesh") != vtkstd::string::npos)
-        {
-        polyMeshFound = true;
-
-        path << "polyMesh/";
-
-        //get number of files in the polyMesh dir
-        vtkDirectory * polyMeshDirectory = vtkDirectory::New();
-        polyMeshDirectory->Open(path.str().c_str());
-        int numPolyMeshFiles = polyMeshDirectory->GetNumberOfFiles();
-        //Look for points/faces files
-        for(int k = 0; k < numPolyMeshFiles; k++)
-          {
-          vtkstd::string tempFile2 = vtkstd::string(polyMeshDirectory->
-                                                    GetFile(k));
-          if(tempFile2.find("points") != vtkstd::string::npos)
-            {
-            this->PolyMeshPointsDir->value[i] = timeStep.str();
-            pointsFound = true;
-            }
-          else if(tempFile2.find("faces") != vtkstd::string::npos)
-            {
-            this->PolyMeshFacesDir->value[i] = timeStep.str();
-            facesFound = true;
-            }
-          }
-
-        //if there is no points or faces found in this timestep
-        //set it equal to previous time step if no previous
-        //set it equal to "constant" dir
-        if(!pointsFound)
-          {
-          if(i != 0)
-            {
-            this->PolyMeshPointsDir->value[i] =
-              this->PolyMeshPointsDir->value[i-1];
-            }
-          else
-            {
-            this->PolyMeshPointsDir->value[i] = vtkstd::string("constant");
-            }
-          }
-        if(!facesFound)
-          {
-          if(i != 0)
-            {
-            this->PolyMeshFacesDir->value[i] =
-              this->PolyMeshFacesDir->value[i-1];
-            }
-          else
-            {
-            this->PolyMeshFacesDir->value[i] = vtkstd::string("constant");
-            }
-          }
-
-        polyMeshDirectory->Delete();
-        break;  //found - stop looking
-        }
+      face& ownerCellFace = cellFaces[ownerCell][--nCellFaces[ownerCell]];
+      ownerCellFace.faceIndex = i;
+      ownerCellFace.neighborFace = false;
       }
 
-    //if there is no polyMesh dir
-    //set  it equal to prev timestep
-    //if no prev set to "constant" dir
-    if(!polyMeshFound)
+    const int neighborCell = faceNeighbor[i]; // must be a signed int
+    if(neighborCell >= 0)
       {
-      if(i != 0)
-        {
-        //set points/faces location to previous timesteps value
-        this->PolyMeshPointsDir->value[i] = this->PolyMeshPointsDir->value[i-1];
-        this->PolyMeshFacesDir->value[i] = this->PolyMeshFacesDir->value[i-1];
-        }
-      else
-        {
-        //set points/faces to constant
-        this->PolyMeshPointsDir->value[i] = vtkstd::string("constant");
-        this->PolyMeshFacesDir->value[i] = vtkstd::string("constant");
-        }
+      face& neighborCellFace
+        = cellFaces[neighborCell][--nCellFaces[neighborCell]];
+      neighborCellFace.faceIndex = i;
+      neighborCellFace.neighborFace = true;
       }
-    directory->Delete();
     }
 
-  vtkDebugMacro(<<"Points/faces list created");
-  return;
+  vtkDebugMacro(<<"Owner file read");
+  return true;
 }
 
 // ****************************************************************************
-//  Method: vtkOpenFOAMReader::GetDataType
+//  Method: vtkOpenFOAMReader::PopulatePolyMeshDirArrays
 //
 //  Purpose:
-//  determines whether a variable is a volume scalar, vector or neither
-//  for meta data
+//  create a Lookup Table containing the location of the points
+//  and faces files for each time steps mesh
 //
 // ****************************************************************************
-const char * vtkOpenFOAMReader::GetDataType(const char * pathIn,
-                                      const char * fileNameIn)
+void vtkOpenFOAMReader::PopulatePolyMeshDirArrays()
 {
-  vtkstd::string path(pathIn);
-  vtkstd::string fileName(fileNameIn);
-  vtkstd::string filePath = path+"/"+fileName;
-  vtkDebugMacro(<<"Get data type of: "<<filePath.c_str());
-  stdString* tempStringStruct;
-  ifstream * input = new ifstream(filePath.c_str(), ios::in VTK_IOS_NOCREATE);
-  //make sure file exists
-  if(input->fail())
-    {
-    input->close();
-    delete input;
-    return "Null";
-    }
-
-  vtkstd::string temp;
-  vtkstd::string foamClass;
-  vtksys_ios::stringstream tokenizer;
-  int opened;
-
-  //see if fileName is a file or directory
-  vtkDirectory * directory = vtkDirectory::New();
-  opened = directory->Open(filePath.c_str());
-  directory->Delete();
-  if(opened)
-    {
-    input->close();
-    delete input;
-    return "Directory";
-    }
+  vtkDebugMacro(<<"Create list of points/faces file directories");
 
-  //find class entry
-  tempStringStruct = this->GetLine(input);
-  temp = tempStringStruct->value;
-  delete tempStringStruct;
-  while(temp.find("class") == vtkstd::string::npos && !input->eof())
-    {
-    tempStringStruct = this->GetLine(input);
-    temp = tempStringStruct->value;
-    delete tempStringStruct;
-    }
+  //intialize size to number of timesteps
+  this->PolyMeshPointsDir->value.resize(this->NumberOfTimeSteps);
+  this->PolyMeshFacesDir->value.resize(this->NumberOfTimeSteps);
 
-  //return type
-  if(!input->eof())
+  //loop through each timestep
+  for(int i = 0; i < this->NumberOfTimeSteps; i++)
     {
-    temp.erase(temp.begin()+temp.find(";"));
-    //PARSE OUT CLASS TYPE
-    tokenizer << temp;
-    //while(tokenizer >> foamClass);
-    while(!tokenizer.eof())
-      {
-      tokenizer >> foamClass;
+    //create the path to the timestep
+    vtkStdString polyMeshPath
+      = *this->PathPrefix + this->TimeNames->GetValue(i) + "/polyMesh/";
+    vtkFoamIOobject io;
+
+    if(io.open(polyMeshPath + "faces") || io.open(polyMeshPath + "faces.gz"))
+      {
+      io.close();
+      //set points/faces location to current timesteps value
+      this->PolyMeshFacesDir->value[i] = this->TimeNames->GetValue(i);
       }
-    //return scalar, vector, or invalid
-    if(foamClass =="volScalarField")
+    else
       {
-      input->close();
-      delete input;
-      return "Scalar";
+      if(i != 0)
+        {
+        //set points/faces location to previous timesteps value
+        this->PolyMeshFacesDir->value[i] = this->PolyMeshFacesDir->value[i-1];
+        }
+      else
+        {
+        //set points/faces to constant
+        this->PolyMeshFacesDir->value[i] = vtkStdString("constant");
+        }
       }
-    else if (foamClass =="volVectorField")
+
+    if(io.open(polyMeshPath + "points") || io.open(polyMeshPath + "points.gz"))
       {
-      input->close();
-      delete input;
-      return "Vector";
+      io.close();
+      //set points/faces location to current timesteps value
+      this->PolyMeshPointsDir->value[i] = this->TimeNames->GetValue(i);
       }
     else
       {
-      input->close();
-      delete input;
-      return "Invalid";
+      if(i != 0)
+        {
+        //set points/faces location to previous timesteps value
+        this->PolyMeshPointsDir->value[i]
+          = this->PolyMeshPointsDir->value[i-1];
+        }
+      else
+        {
+        //set points/faces to constant
+        this->PolyMeshPointsDir->value[i] = vtkStdString("constant");
+        }
       }
     }
-
-  //if the file format is wrong return invalid
-  else
-    {
-    input->close();
-    delete input;
-    return "invalid";
-    }
+  vtkDebugMacro(<<"Points/faces list created");
+  return;
 }
 
 // ****************************************************************************
-//  Method: vtkOpenFOAMReader::GetInternalVariableAtTimestep
+//  Method: vtkOpenFOAMReader::GetVariableAtTimestep
 //
 //  Purpose:
-//  returns the values for a request variable for the internal mesh
+//  returns the values for a request variable for the internal and boundary
+//  meshes
 //
 // ****************************************************************************
-vtkDoubleArray * vtkOpenFOAMReader::GetInternalVariableAtTimestep
-                 (const char * varNameIn, int timeState)
+bool vtkOpenFOAMReader::GetVariableAtTimestep(
+  vtkUnstructuredGrid* internalMesh, unstructuredGridVector* boundaryMesh,
+  vtkFoamDict* boundaryDictPtr, const char* varNameIn, int timeState)
 {
-  vtkstd::string varName(varNameIn);
-  vtksys_ios::stringstream varPath;
-  varPath << this->PathPrefix->value << this->Steps[timeState] << "/" << varName;
-  vtkDebugMacro(<<"Get internal variable: "<<varPath.str().c_str());
-  vtkDoubleArray *data = vtkDoubleArray::New();
+  vtkStdString varPath = *this->PathPrefix
+    + this->TimeNames->GetValue(timeState) + "/" + vtkStdString(varNameIn);
+  vtkDebugMacro(<<"Get variable: "<<varPath.c_str());
 
-  vtkstd::string temp;
-  stdString* tempStringStruct;
-  bool binaryWriteFormat;
-  ifstream * input = new ifstream(varPath.str().c_str(), ios::in VTK_IOS_NOCREATE);
-  //make sure file exists
-  if(input->fail())
+  // open the file
+  vtkFoamIOobject io;
+  if(!io.open(varPath))
     {
-    input->close();
-    delete input;
-    return data;
+    vtkErrorMacro(<<"Error opening " << io.fileName().c_str() << ": "
+      << io.error().str().c_str());
+    return false;
     }
 
-  //determine if file is binary or ascii
-  while(temp.find("format") == vtkstd::string::npos)
+  // if the variable is disabled on selection panel then skip it
+  if(this->CellDataArraySelection->ArrayExists(io.objectName().c_str())
+    && !this->CellDataArraySelection->ArrayIsEnabled(io.objectName().c_str()))
     {
-    tempStringStruct = this->GetLine(input);
-    temp = tempStringStruct->value;
-    delete tempStringStruct;
+    return true;
     }
-  input->close();
 
-  //reopen file in correct format
-  if(temp.find("binary") != vtkstd::string::npos)
-    {
-#ifdef _WIN32
-    input->open(varPath.str().c_str(), ios::binary | ios::in VTK_IOS_NOCREATE);
-#else
-    input->open(varPath.str().c_str(), ios::in VTK_IOS_NOCREATE);
-#endif
-    binaryWriteFormat = true;
-    }
-  else
+  // read the field file into dictionary
+  vtkFoamDict dict;
+  if(!dict.read(io))
     {
-    input->open(varPath.str().c_str(),ios::in);
-    binaryWriteFormat = false;
+    vtkErrorMacro(<<"Error reading line " << io.lineNumber()
+      << " of " << io.fileName().c_str() << ": " << io.error().str().c_str());
+    return false;
     }
 
-  vtkstd::string foamClass;
-  vtksys_ios::stringstream tokenizer;
-  double value;
-
-  //find class
-  tempStringStruct = this->GetLine(input);
-  temp = tempStringStruct->value;
-  delete tempStringStruct;
-
-  while(temp.find("class") == vtkstd::string::npos)
+  if(dict.type() != vtkFoamDict::DICTIONARY)
     {
-    tempStringStruct = this->GetLine(input);
-    temp = tempStringStruct->value;
-    delete tempStringStruct;
+    vtkErrorMacro(<<"File " << io.fileName().c_str()
+      << "is not valid as a field file");
+    return false;
     }
-  temp.erase(temp.begin()+temp.find(";"));
-  tokenizer << temp;
-  //while(tokenizer >> foamClass);
-  while(!tokenizer.eof())
+
+  // set internal values
+  vtkDoubleArray* iData;
+  vtkFoamEntry& iEntry = dict.lookup("internalField");
+  if(!iEntry.found())
     {
-    tokenizer >> foamClass;
+    vtkErrorMacro(<<"internalField not found in " << io.fileName().c_str());
+    return false;
     }
-  temp="";
-  //create scalar arrays
-  if(foamClass =="volScalarField")
+  if(iEntry.firstValue().type() == vtkFoamEntryValue::EMPTYLIST)
     {
-    while(temp.find("internalField") == vtkstd::string::npos)
+    if(this->NumCells == 0)
       {
-      tempStringStruct = this->GetLine(input);
-      temp = tempStringStruct->value;
-      delete tempStringStruct;
+      // if there's no cell there shouldn't be any boundary faces either
+      return true;
       }
-    //nonuniform
-    if(!(temp.find("nonuniform") == vtkstd::string::npos))
+    else
       {
-      //create an array
-      tempStringStruct = this->GetLine(input);
-      temp = tempStringStruct->value;
-      delete tempStringStruct;
-
-      int scalarCount;
-      tokenizer.clear();
-      tokenizer << temp;
-      tokenizer >> scalarCount;
-      data->SetNumberOfValues(NumCells);
-
-      //binary data
-      if(binaryWriteFormat)
-        {
-        //add values to array
-        input->get(); //parenthesis
-        for(int i = 0; i < scalarCount; i++)
-          {
-          input->read((char *) &value, sizeof(double));
-          data->SetValue(i, value);
-          }
-        }
+      vtkErrorMacro(<<"internalField is empty");
+      return false;
+      }
+    }
 
-      //ascii data
-      else
+  if(io.className() == "volScalarField")
+    {
+    if(iEntry.firstValue().isUniform())
+      {
+      iData = vtkDoubleArray::New();
+      iData->SetNumberOfValues(this->NumCells);
+      double num;
+      iEntry >> num;
+      for(int i = 0; i < this->NumCells; i++)
         {
-        //add values to array
-        tempStringStruct = this->GetLine(input);
-        temp = tempStringStruct->value;
-        delete tempStringStruct; //discard (
-
-        for(int i = 0; i < scalarCount; i++)
-          {
-          tempStringStruct = this->GetLine(input);
-          temp = tempStringStruct->value;
-          delete tempStringStruct;
-          tokenizer.clear();
-          tokenizer << temp;
-          tokenizer >> value;
-          data->SetValue(i, value);
-          }
+        iData->SetValue(i, num);
         }
       }
-
-    //uniform
-    else if(!(temp.find("uniform") == vtkstd::string::npos))
+    else // nonuniform
       {
-      //parse out the uniform value
-      vtkstd::string token;
-      temp.erase(temp.begin()+temp.find(";"));
-      tokenizer.clear();
-      tokenizer << temp;
-      //while(tokenizer>>token);
-      while(!tokenizer.eof())
+      if(iEntry.firstValue().type() != vtkFoamEntryValue::SCALARLIST)
         {
-        tokenizer >> token;
+        vtkErrorMacro(<<"internalField is not a scalarField");
+        return false;
         }
-      tokenizer.clear();
-      tokenizer << token;
-      tokenizer >> value;
-      data->SetNumberOfValues(NumCells);
-
-      //create array of uniform values
-      for(int i = 0; i < NumCells; i++)
+      if(iEntry.scalarList().GetSize() != this->NumCells)
         {
-        data->SetValue(i, value);
+        vtkErrorMacro(<<"Number of cells in mesh and field don't match: "
+          << "mesh = " << this->NumCells << ", field = "
+          << iEntry.scalarList().GetSize());
+        return false;
         }
-      }
-
-    //no data
-    else
-      {
-      input->close();
-      delete input;
-      return data;
+      iData = (vtkDoubleArray *)iEntry.ptr();
       }
     }
-
-  //create vector arrays
-  else if(foamClass == "volVectorField")
+  else if(io.className() == "volVectorField")
     {
-    tempStringStruct = this->GetLine(input);
-    temp = tempStringStruct->value;
-    delete tempStringStruct;
-    while(temp.find("internalField") == vtkstd::string::npos)
+    if(iEntry.firstValue().isUniform())
       {
-      tempStringStruct = this->GetLine(input);
-      temp = tempStringStruct->value;
-      delete tempStringStruct;
-      }
-    if(!(temp.find("nonuniform") == vtkstd::string::npos))
-      {
-      //create an array
-      tempStringStruct = this->GetLine(input);
-      temp = tempStringStruct->value;
-      delete tempStringStruct;
-
-      int vectorCount;
-      tokenizer.clear();
-      tokenizer << temp;
-      tokenizer >> vectorCount;
-      data->SetNumberOfComponents(3);
-
-      //binary data
-      if(binaryWriteFormat)
+      iData = vtkDoubleArray::New();
+      iData->SetNumberOfComponents(3);
+      iData->SetNumberOfTuples(this->NumCells);
+      // have to determine the type of vector
+      if(iEntry.firstValue().type() == vtkFoamEntryValue::LABELLIST)
         {
-        //add values to the array
-        input->get(); //parenthesis
-        for(int i = 0; i < vectorCount; i++)
+        const vtkstd::vector<int>& ll = iEntry.labelList();
+        if(ll.size() != 3)
           {
-          input->read((char *) &value, sizeof(double));
-          data->InsertComponent(i, 0, value);
-          input->read((char *) &value, sizeof(double));
-          data->InsertComponent(i, 1, value);
-          input->read((char *) &value, sizeof(double));
-          data->InsertComponent(i, 2, value);
+          vtkErrorMacro(<<"The uniform list isn't a vector of size 3");
+          iData->Delete();
+          return false;
+          }
+        double v[3];
+        v[0] = static_cast<double>(ll[0]);
+        v[1] = static_cast<double>(ll[1]);
+        v[2] = static_cast<double>(ll[2]);
+        for(int i = 0; i < this->NumCells; i++)
+          {
+          iData->SetTuple(i, v);
           }
         }
-
-      //ascii data
-      else
+      else if(iEntry.firstValue().type() == vtkFoamEntryValue::SCALARLIST)
         {
-        //add values to the array
-        tempStringStruct = this->GetLine(input);
-        temp = tempStringStruct->value;
-        delete tempStringStruct; //discard (
-        for(int i = 0; i < vectorCount; i++)
+        vtkDoubleArray& sl = iEntry.scalarList();
+        if(sl.GetSize() != 3)
           {
-          tempStringStruct = this->GetLine(input);
-          temp = tempStringStruct->value;
-          delete tempStringStruct;
-
-          //REMOVE BRACKETS
-          temp.erase(temp.begin()+temp.find("("));
-          temp.erase(temp.begin()+temp.find(")"));
-
-          //GRAB X,Y,&Z VALUES
-          tokenizer.clear();
-          tokenizer << temp;
-          tokenizer >> value;
-          data->InsertComponent(i, 0, value);
-          tokenizer >> value;
-          data->InsertComponent(i, 1, value);
-          tokenizer >> value;
-          data->InsertComponent(i, 2, value);
+          vtkErrorMacro(
+            << "The uniform list isn't a vector of size 3 (actual size "
+            << sl.GetSize() << ")");
+          iData->Delete();
+          return false;
+          }
+        double* v = sl.GetPointer(0);
+        for(int i = 0; i < this->NumCells; i++)
+          {
+          iData->SetTuple(i, v);
           }
         }
-      }
-    else if(!(temp.find("uniform") == vtkstd::string::npos))
-      {
-      //create an array of uniform values
-      double value1, value2, value3;
-
-      //parse out the uniform values
-      temp.erase(temp.begin(), temp.begin()+temp.find("(")+1);
-      temp.erase(temp.begin()+temp.find(")"), temp.end());
-      tokenizer.clear();
-      tokenizer << temp;
-      tokenizer >> value1;
-      tokenizer >> value2;
-      tokenizer >> value3;
-      data->SetNumberOfComponents(3);
-      for(int i = 0; i < NumCells; i++)
+      else
         {
-        data->InsertComponent(i, 0, value1);
-        data->InsertComponent(i, 1, value2);
-        data->InsertComponent(i, 2, value3);
+        vtkErrorMacro(<<"Wrong list type");
+        iData->Delete();
+        return false;
         }
       }
-
-    //no data
-    else
+    else // nonuniform
       {
-      input->close();
-      delete input;
-      return data;
+      if(iEntry.firstValue().type() != vtkFoamEntryValue::VECTORLIST)
+        {
+        vtkErrorMacro(<<"internalField is not a vectorField");
+        return false;
+        }        
+      if(iEntry.vectorList().GetNumberOfTuples() != this->NumCells)
+        {
+        vtkErrorMacro(<<"Number of cells in mesh and field don't match: "
+          << "mesh = " << this->NumCells << ", field = "
+          << iEntry.vectorList().GetNumberOfTuples());
+        return false;
+        }
+      iData = (vtkDoubleArray *)iEntry.ptr();
       }
     }
-  input->close();
-  delete input;
-  vtkDebugMacro(<<"Internal variable data read");
-  return data;
-}
-
-// ****************************************************************************
-//  Method: vtkOpenFOAMReader::GetBoundaryVariableAtTimestep
-//
-//  Purpose:
-//  returns the values for a request variable for a bondary region
-//
-// ****************************************************************************
-vtkDoubleArray * vtkOpenFOAMReader::GetBoundaryVariableAtTimestep
-                (int boundaryIndex, const char * varNameIn, int timeState,
-                 vtkUnstructuredGrid * internalMesh)
-{
-  vtkstd::string varName(varNameIn);
-  vtksys_ios::stringstream varPath;
-  varPath << this->PathPrefix->value << this->Steps[timeState] << "/" << varName;
-  vtkDebugMacro(<<"Get boundary variable: "<<varPath.str().c_str());
-  vtkDoubleArray *data = vtkDoubleArray::New();
-
-  vtkstd::string temp;
-  stdString* tempStringStruct;
-  bool binaryWriteFormat;
-  ifstream * input = new ifstream(varPath.str().c_str(), ios::in VTK_IOS_NOCREATE);
-  //make sure file exists
-  if(input->fail())
+  else
     {
-    input->close();
-    delete input;
-    return data;
+    vtkErrorMacro(<<"Unsupported field class");
+    return false;
     }
 
-  //determine if file is binary or ascii
-  while(temp.find("format") == vtkstd::string::npos)
+  // Add field only if internal Mesh exists (skip if not selected)
+  // Note we still need to read internalField even if internal mesh is
+  // not selected, since boundaries without value entries may refer to
+  // the internalField.
+  if((iData->GetSize() > 0) 
+    && (this->GetPatchArrayStatus(this->GetPatchArrayName(0)))
+    && (internalMesh != NULL))
     {
-    tempStringStruct = this->GetLine(input);
-    temp = tempStringStruct->value;
-    delete tempStringStruct;
+    iData->SetName(io.objectName().c_str());
+    internalMesh->GetCellData()->AddArray(iData);
     }
-  input->close();
 
-  //reopen file in correct format
-  if(temp.find("binary") != vtkstd::string::npos)
-    {
-#ifdef _WIN32
-    input->open(varPath.str().c_str(), ios::binary | ios::in VTK_IOS_NOCREATE);
-#else
-    input->open(varPath.str().c_str(), ios::in VTK_IOS_NOCREATE);
-#endif
-    binaryWriteFormat = true;
-    }
-  else
+  // set boundary values
+  vtkFoamEntry& bEntry = dict.lookup("boundaryField");
+  vtkFoamDict& boundaryDict = *boundaryDictPtr;
+  for(size_t boundaryIndex = 0, i = 0; i < boundaryDict.size(); i++)
     {
-    input->open(varPath.str().c_str(),ios::in);
-    binaryWriteFormat = false;
-    }
-
-  vtkstd::string foamClass;
-  vtksys_ios::stringstream tokenizer;
-  double value;
+    vtkFoamEntry& boundaryEntryI = boundaryDict.entry(i);
+    int nFaces;
+    boundaryEntryI.dictionary().lookup("nFaces") >> nFaces;
 
-  //find class
-  tempStringStruct = this->GetLine(input);
-  temp = tempStringStruct->value;
-  delete tempStringStruct;
-  while(temp.find("class") == vtkstd::string::npos)
-    {
-    tempStringStruct = this->GetLine(input);
-    temp = tempStringStruct->value;
-    delete tempStringStruct;
-    }
-  temp.erase(temp.begin()+temp.find(";"));
-  tokenizer << temp;
-  //while(tokenizer >> foamClass);
-  while(!tokenizer.eof())
-    {
-    tokenizer >> foamClass;
-    }
-  temp="";
-  //create scalar arrays
-  if(foamClass =="volScalarField")
-    {
-    //find desired mesh
-    while(temp.find(this->BoundaryNames->value[boundaryIndex]) == 
-            vtkstd::string::npos && !input->eof())
+    // skip empty boundary
+    if(nFaces == 0)
       {
-      tempStringStruct = this->GetLine(input);
-      temp = tempStringStruct->value;
-      delete tempStringStruct;
+      continue;
       }
-    if(input->eof())
+
+    const vtkStdString& boundaryNameI = boundaryEntryI.keyword();
+
+    // Add field only if given boundary is selected for display
+    if(!this->GetPatchArrayStatus(boundaryNameI.c_str()))
       {
-      input->close();
-      delete input;
-      return data;
+      continue;
       }
-    //find value entry
-    while(temp.find("}") == vtkstd::string::npos &&
-          temp.find("value ") == vtkstd::string::npos)
+
+    vtkFoamEntry& bEntryI = bEntry.dictionary().lookup(boundaryNameI);
+    if(!bEntryI.found())
       {
-      tempStringStruct = this->GetLine(input);
-      temp = tempStringStruct->value;
-      delete tempStringStruct; //find value
+      vtkErrorMacro(<< "boundaryField " << boundaryNameI.c_str()
+        << " not found in object " << varNameIn << " at time = "
+        << this->TimeNames->GetValue(timeState).c_str());
+      iData->Delete();
+      return false;
       }
 
-    //nonuniform
-    if(!(temp.find("nonuniform") == vtkstd::string::npos))
+    if(bEntryI.firstValue().type() != vtkFoamEntryValue::DICTIONARY)
       {
+      vtkErrorMacro(<< "Type of boundaryField " << boundaryNameI.c_str()
+        << " is not a subdictionary in object " << varNameIn << " at time = "
+        << this->TimeNames->GetValue(timeState).c_str());
+      iData->Delete();
+      return false;
+      }
 
-      //binary data
-      if(binaryWriteFormat)
-        {
-        //create an array
-        tempStringStruct = this->GetLine(input);
-        temp = tempStringStruct->value;
-        delete tempStringStruct;
-
-        int scalarCount;
-        tokenizer.clear();
-        tokenizer << temp;
-        tokenizer >> scalarCount;
-        data->SetNumberOfValues(scalarCount);
+    vtkDoubleArray* vData;
 
-        //assign values to the array
-        input->get(); //parenthesis
-        for(int i = 0; i < scalarCount; i++)
+    vtkFoamEntry& vEntry = bEntryI.dictionary().lookup("value");
+    // uniformValue is for Hrv's dev version
+    vtkFoamEntry& uvEntry = bEntryI.dictionary().lookup("uniformValue");
+    
+    if(vEntry.found() || uvEntry.found()) // the boundary has a value entry
+      {
+      if(io.className() == "volScalarField")
+        {
+        if((vEntry.found() && vEntry.firstValue().isUniform())
+          || uvEntry.found())
+          {
+          vData = vtkDoubleArray::New();
+          vData->SetNumberOfValues(nFaces);
+          double num;
+          if(vEntry.found()) // "value" entry has the priority
+            {
+            vEntry >> num;
+            }
+          else
+            {
+            uvEntry >> num;
+            }
+          for(int j = 0; j < nFaces; j++)
+            {
+            vData->SetValue(j, num);
+            }
+          }
+        else // nonuniform
           {
-          input->read((char *) &value, sizeof(double));
-          data->SetValue(i, value);
+          if(vEntry.firstValue().type() != vtkFoamEntryValue::SCALARLIST)
+            {
+            vtkErrorMacro(<<"boundaryField " << boundaryNameI.c_str()
+              << " is not a scalarField");
+            iData->Delete();
+            return false;
+            }
+          if(vEntry.scalarList().GetSize() != nFaces)
+            {
+            vtkErrorMacro(<<"Number of faces in field and bondary don't match");
+            iData->Delete();
+            return false;
+            }
+          vData = (vtkDoubleArray *)vEntry.ptr();
           }
         }
-
-      //ascii data
-      else
+      else if(io.className() == "volVectorField")
         {
-        temp.erase(temp.begin(), temp.begin()+temp.find(">")+1);
-        //ascii data with 10 or less values are on the same line
-        //>10
-        if(temp == vtkstd::string(" "))
-          {
-          //create an array of data
-          tempStringStruct = this->GetLine(input);
-          temp = tempStringStruct->value;
-          delete tempStringStruct;
-
-          int scalarCount;
-          tokenizer.clear();
-          tokenizer << temp;
-          tokenizer >> scalarCount;
-          data->SetNumberOfValues(scalarCount);
-          tempStringStruct = this->GetLine(input);
-          temp = tempStringStruct->value;
-          delete tempStringStruct; //discard (
-
-          for(int i = 0; i < scalarCount; i++)
-            {
-            tempStringStruct = this->GetLine(input);
-            temp = tempStringStruct->value;
-            delete tempStringStruct;
-
-            tokenizer.clear();
-            tokenizer << temp;
-            tokenizer >> value;
-            data->SetValue(i, value);
+        if((vEntry.found() && vEntry.firstValue().isUniform())
+          || uvEntry.found())
+          {
+          vData = vtkDoubleArray::New();
+          vData->SetNumberOfComponents(3);
+          vData->SetNumberOfTuples(nFaces);
+          // the "value" entry has the priority
+          vtkFoamEntry& vvEntry = vEntry.found() ? vEntry : uvEntry;
+          // have to determine the type of vector
+          if(vvEntry.firstValue().type() == vtkFoamEntryValue::LABELLIST)
+            {
+            const vtkstd::vector<int>& ll = vvEntry.labelList();
+            if(ll.size() != 3)
+              {
+              vtkErrorMacro(<<"The uniform list isn't a vector of size 3");
+              iData->Delete();
+              vData->Delete();
+              return false;
+              }
+            double v[3];
+            v[0] = static_cast<double>(ll[0]);
+            v[1] = static_cast<double>(ll[1]);
+            v[2] = static_cast<double>(ll[2]);
+            for(int j = 0; j < nFaces; j++)
+              {
+              vData->SetTuple(j, v);
+              }
+            }
+          else if(vvEntry.firstValue().type() == vtkFoamEntryValue::SCALARLIST)
+            {
+            vtkDoubleArray& sl = vvEntry.scalarList();
+            if(sl.GetSize() != 3)
+              {
+              vtkErrorMacro(<< "The uniform value list in patch "
+                << boundaryNameI.c_str() << " of file "
+                << io.fileName().c_str()
+                << " isn't a vector of size 3 (actual size " << sl.GetSize()
+                << ")");
+              iData->Delete();
+              vData->Delete();
+              return false;
+              }
+            double* v = sl.GetPointer(0);
+            for(int j = 0; j < nFaces; j++)
+              {
+              vData->SetTuple(j, v);
+              }
+            }
+          else
+            {
+            vtkErrorMacro(<<"Wrong list type");
+            iData->Delete();
+            vData->Delete();
+            return false;
             }
           }
-        //=<10
-        else
+        else // nonuniform
           {
-          //create an array with data
-          int scalarCount;
-          tokenizer.clear();
-          tokenizer << temp;
-          tokenizer >> scalarCount;
-          data->SetNumberOfValues(scalarCount);
-          temp.erase(temp.begin(), temp.begin()+temp.find("(")+1);
-          temp.erase(temp.begin()+temp.find(")"), temp.end());
-
-          tokenizer.clear();
-          tokenizer << temp;
-          for(int i = 0; i < scalarCount; i++)
+          if(vEntry.firstValue().type() != vtkFoamEntryValue::VECTORLIST)
+            {
+            vtkErrorMacro(<<"boundaryField " << boundaryNameI.c_str()
+              << " is not a vectorField");
+            iData->Delete();
+            return false;
+            }
+          if(vEntry.vectorList().GetNumberOfTuples() != nFaces)
             {
-            tokenizer >> value;
-            data->SetValue(i, value);
+            vtkErrorMacro(<<"Number of faces in field and bondary don't match");
+            iData->Delete();
+            return false;
             }
+          vData = (vtkDoubleArray *)vEntry.ptr();
           }
         }
-      }
-
-    //uniform
-    else if(!(temp.find("uniform") == vtkstd::string::npos))
-      {
-      //create an array of uniform values
-      double value1 = 0;
-      temp.erase(temp.begin(), temp.begin()+temp.find("uniform")+7);
-      temp.erase(temp.begin()+temp.find(";"), temp.end());
-      tokenizer.clear();
-      tokenizer << temp;
-      tokenizer >> value1;
-      data->SetNumberOfValues(this->SizeOfBoundary->value[boundaryIndex]);
-      for(int i = 0; i < this->SizeOfBoundary->value[boundaryIndex]; i++)
-        {
-        data->SetValue(i, value1);
-        }
-      }
-
-    //no data
-    else
-      {
-      int cellId;
-      vtkDataArray * internalData = internalMesh->GetCellData()->
-                                    GetArray(varName.c_str());
-      data->SetNumberOfValues(this->SizeOfBoundary->value[boundaryIndex]);
-      for(int i = 0; i < this->SizeOfBoundary->value[boundaryIndex]; i++)
+      else // other field types
         {
-        cellId = this->FaceOwner->GetValue(this->StartFace + i);
-        data->SetValue(i, internalData->GetComponent(cellId, 0));
+        vtkErrorMacro(<<"Unsupported field class");
+        iData->Delete();
+        return false;
         }
-      input->close();
-      delete input;
-      return data;
-      }
-    }
-  //CREATE VECTOR ARRAYS
-  else if(foamClass == "volVectorField")
-    {
-    while(temp.find(this->BoundaryNames->value[boundaryIndex]) == 
-          vtkstd::string::npos && !input->eof())
-      {
-      tempStringStruct = this->GetLine(input);
-      temp = tempStringStruct->value;
-      delete tempStringStruct;
-      }
-    if(input->eof())
-      {
-      input->close();
-      delete input;
-      return data;
       }
-    while(temp.find("}") == vtkstd::string::npos &&
-          temp.find("value ") == vtkstd::string::npos)
+    else // doesn't have a value entry
       {
-      tempStringStruct = this->GetLine(input);
-      temp = tempStringStruct->value;
-      delete tempStringStruct; //find value
-      }
-    //nonuniform
-    if(!(temp.find("nonuniform") == vtkstd::string::npos))
-      {
-      //create an array
-      tempStringStruct = this->GetLine(input);
-      temp = tempStringStruct->value;
-      delete tempStringStruct;
+      int startFace;
+      boundaryEntryI.dictionary().lookup("startFace") >> startFace;
+      vData = vtkDoubleArray::New();
 
-      int vectorCount;
-      tokenizer.clear();
-      tokenizer << temp;
-      tokenizer >> vectorCount;
-      data->SetNumberOfComponents(3);
-
-      //binary data
-      if(binaryWriteFormat)
+      if(io.className() == "volScalarField")
         {
-        //insert values into the array
-        input->get(); //parenthesis
-        for(int i = 0; i < vectorCount; i++)
+        vData->SetNumberOfValues(nFaces);
+        for(int j = 0; j < nFaces; j++)
           {
-          input->read((char *) &value, sizeof(double));
-          data->InsertComponent(i, 0, value);
-          input->read((char *) &value, sizeof(double));
-          data->InsertComponent(i, 1, value);
-          input->read((char *) &value, sizeof(double));
-          data->InsertComponent(i, 2, value);
+          const int cellId = this->FaceOwner->GetValue(startFace + j);
+          vData->SetValue(j, iData->GetValue(cellId));
           }
         }
-
-        //ascii data
-        else
+      else if(io.className() == "volVectorField")
+        {
+        vData->SetNumberOfComponents(3);
+        vData->SetNumberOfTuples(nFaces);
+        for(int j = 0; j < nFaces; j++)
           {
-          //insert values into the array
-          tempStringStruct = this->GetLine(input);
-          temp = tempStringStruct->value;
-          delete tempStringStruct; //discard (
-          for(int i = 0; i < vectorCount; i++)
-            {
-            tempStringStruct = this->GetLine(input);
-            temp = tempStringStruct->value;
-            delete tempStringStruct;
-
-            //REMOVE BRACKETS
-            temp.erase(temp.begin()+temp.find("("));
-            temp.erase(temp.begin()+temp.find(")"));
-
-            //GRAB X,Y,&Z VALUES
-            tokenizer.clear();
-            tokenizer << temp;
-            tokenizer >> value;
-            data->InsertComponent(i, 0, value);
-            tokenizer >> value;
-            data->InsertComponent(i, 1, value);
-            tokenizer >> value;
-            data->InsertComponent(i, 2, value);
-            }
+          const int cellId = this->FaceOwner->GetValue(startFace + j);
+          // in this case GetTuple should be thread-safe accoding to
+          // the FAQ on ParaView Wiki
+          vData->SetTuple(j, iData->GetTuple(cellId));
           }
         }
-
-    //uniform
-    else if(!(temp.find("uniform") == vtkstd::string::npos))
-      {
-      //create an array of uniform values
-      double value1 = 0, value2 = 0, value3 = 0;
-      temp.erase(temp.begin(), temp.begin()+temp.find("(")+1);
-      temp.erase(temp.begin()+temp.find(")"), temp.end());
-      tokenizer.clear();
-      tokenizer << temp;
-      tokenizer >> value1;
-      tokenizer >> value2;
-      tokenizer >> value3;
-
-      data->SetNumberOfComponents(3);
-      for(int i = 0; i < this->SizeOfBoundary->value[boundaryIndex]; i++)
+      else
         {
-        data->InsertComponent(i, 0, value1);
-        data->InsertComponent(i, 1, value2);
-        data->InsertComponent(i, 2, value3);
+        vtkErrorMacro(<<"Unsupported field class");
+        iData->Delete();
+        vData->Delete();
+        return false;
         }
       }
-
-    //no data
-    else
+    if(vData->GetSize() > 0)
       {
-      int cellId;
-      vtkDataArray * internalData = internalMesh->GetCellData()->
-                                    GetArray(varName.c_str());
-      data->SetNumberOfComponents(3);
-      for(int i = 0; i < this->SizeOfBoundary->value[boundaryIndex]; i++)
-        {
-        cellId = this->FaceOwner->GetValue(this->StartFace + i);
-        data->InsertComponent(i, 0, internalData->GetComponent(cellId, 0));
-        data->InsertComponent(i, 1, internalData->GetComponent(cellId, 1));
-        data->InsertComponent(i, 2, internalData->GetComponent(cellId, 2));
-        }
-      input->close();
-      delete input;
-      return data;
+      vData->SetName(io.objectName().c_str());
+      boundaryMesh->value[boundaryIndex]->GetCellData()
+        ->AddArray(vData);
       }
+    vData->Delete();
+    boundaryIndex++;
     }
-  input->close();
-  delete input;
-  vtkDebugMacro(<<"Boundary data read");
-  return data;
+
+  iData->Delete();
+  vtkDebugMacro(<<"Internal variable data read");
+  return true;
 }
 
 // ****************************************************************************
 //  Method: vtkOpenFOAMReader::GatherBlocks
 //
 //  Purpose:
-//  returns a vector of block names for a specified domain
+//  returns a dictionary of block names for a specified domain
 //
 // ****************************************************************************
-stringVector * vtkOpenFOAMReader::GatherBlocks(const char * typeIn,
-                                              int timeState)
+vtkOpenFOAMReader::vtkFoamDict* vtkOpenFOAMReader::GatherBlocks(
+  const char* typeIn, int timeState, bool mustRead)
 {
-  vtkstd::string type(typeIn);
-  vtkstd::string blockPath = this->PathPrefix->value +
-                             this->PolyMeshFacesDir->value[timeState] +
-                             "/polyMesh/"+type;
-  vtkstd::vector< vtkstd::string > blocks;
-  stringVector *returnValue = new stringVector;
-  vtkDebugMacro(<<"Get blocks: "<<blockPath.c_str());
-
-  ifstream * input = new ifstream(blockPath.c_str(), ios::in VTK_IOS_NOCREATE);
-  //if there is no file return a null vector
-  if(input->fail())
-    {
-    input->close();
-    delete input;
-    returnValue->value = blocks;
-    return returnValue;
-    }
-
-  vtkstd::string temp;
-  stdString* tempStringStruct;
-  vtkstd::string token;
-  vtksys_ios::stringstream tokenizer;
-  vtkstd::string tempName;
-
-  //find end of header
-  //while(temp.compare(0,4,"// *",0,4)!=0)
-  while (strcmp(temp.substr(0,4).c_str(), "// *"))
-    {
-    tempStringStruct = this->GetLine(input);
-    temp = tempStringStruct->value;
-    delete tempStringStruct;
-    }
-  tempStringStruct = this->GetLine(input);
-  temp = tempStringStruct->value;
-  delete tempStringStruct;
-
-  tempStringStruct = this->GetLine(input);
-  temp = tempStringStruct->value;
-  delete tempStringStruct; //throw out blank line
-
-  //Number of blocks
-  tokenizer << temp;
-  tokenizer >> this->NumBlocks;
-  blocks.resize(this->NumBlocks);
-
-  //loop through each block
-  for(int i = 0; i < this->NumBlocks; i++)
-    {
-    tempStringStruct = this->GetLine(input);
-    temp = tempStringStruct->value;
-    delete tempStringStruct; //throw out blank line
-
-    //NAME
-    tempStringStruct = this->GetLine(input);
-    temp = tempStringStruct->value;
-    delete tempStringStruct; //name
-
-    tokenizer.clear();
-    tokenizer << temp;
-    tokenizer >> tempName;
-    blocks[i] = tempName;
-    //while(temp.compare(0,1,"}",0,1) != 0)
-    while (strcmp(temp.substr(0,1).c_str(), "}"))
-      {
-      tempStringStruct = this->GetLine(input);
-      temp = tempStringStruct->value;
-      delete tempStringStruct;
-      }
-    }
-  returnValue->value = blocks;
-  input->close();
-  delete input;
-  return returnValue;
+  vtkStdString type(typeIn);
+  vtkStdString blockPath = *this->PathPrefix
+    + this->PolyMeshFacesDir->value[timeState] + "/polyMesh/" + type;
+
+  vtkFoamIOobject io;
+  if(!(io.open(blockPath) || io.open(blockPath + ".gz")))
+    {
+    if(mustRead)
+      {
+      vtkErrorMacro(<<"Error opening " << io.fileName().c_str() << ": "
+        << io.error().str().c_str());
+      }
+    return NULL;
+    }
+
+  vtkFoamDict* dictPtr = new vtkFoamDict;
+  vtkFoamDict& dict = *dictPtr;
+  if(!dict.read(io))
+    {
+    vtkErrorMacro(<<"Error reading line " << io.lineNumber()
+      << " of " << io.fileName().c_str() << ": " << io.error().str().c_str());
+    delete dictPtr;
+    return NULL;
+    }
+  if(dict.type() != vtkFoamDict::DICTIONARY)
+    {
+    vtkErrorMacro(<<"The file type of " << io.fileName().c_str()
+      << " is not a dictionary");
+    delete dictPtr;
+    return NULL;
+    }
+  return dictPtr;
 }
 
 // ****************************************************************************
@@ -2378,138 +4350,105 @@
 //  returns a requested boundary mesh
 //
 // ****************************************************************************
-vtkUnstructuredGrid * vtkOpenFOAMReader::GetBoundaryMesh(int timeState,
-                                                         int boundaryIndex)
+unstructuredGridVector* vtkOpenFOAMReader::GetBoundaryMesh(
+  vtkFoamDict* boundaryDictPtr, intVectorVector* facesPoints,
+  vtkPoints* points)
 {
-  vtkUnstructuredGrid * boundaryMesh = vtkUnstructuredGrid::New();
-  vtkstd::string boundaryPath = this->PathPrefix->value +
-                                this->PolyMeshFacesDir->value[timeState] +
-                                "/polyMesh/boundary";
-  vtkDebugMacro(<<"Create boundary mesh: "<<boundaryPath.c_str());
+  vtkFoamDict& boundaryDict = *boundaryDictPtr;
+  unstructuredGridVector* boundaryMesh = new unstructuredGridVector;
 
-  int nFaces;
-
-  ifstream * input = new ifstream(boundaryPath.c_str(), ios::in VTK_IOS_NOCREATE);
-  //return a Null object
-  if(input->fail())
+  for(size_t boundaryIndex = 0; boundaryIndex < boundaryDict.size();
+    boundaryIndex++)
     {
-    input->close();
-    delete input;
-    return boundaryMesh;
-    }
+    vtkFoamDict& dict = boundaryDict.entry(boundaryIndex).dictionary();
 
-  vtkstd::string temp;
-  stdString* tempStringStruct;
-  vtkstd::string token;
-  vtksys_ios::stringstream tokenizer;
+    int nFaces;
+    vtkFoamEntry& nFacesEntry = dict.lookup("nFaces");
+    vtkFoamEntry& startFaceEntry = dict.lookup("startFace");
+    if(!nFacesEntry.found() || !startFaceEntry.found())
+      {
+      vtkWarningMacro(
+        << "Entry nFaces or startFace not found in boundary entry "
+        << boundaryDict.entry(boundaryIndex).keyword().c_str());
+      for(size_t i = 0; i < boundaryMesh->value.size(); i++)
+        {
+        boundaryMesh->value[i]->Delete();
+        }
+      delete boundaryMesh;
+      return NULL;
+      }
+    nFacesEntry >> nFaces;
 
-  //find desired mesh entry
-  while(temp.find(this->BoundaryNames->value[boundaryIndex]) == 
-        vtkstd::string::npos)
-    {
-    tempStringStruct = this->GetLine(input);
-    temp = tempStringStruct->value;
-    delete tempStringStruct;
-    }
+    // skip empty boundary
+    if(nFaces == 0)
+      {
+      continue;
+      }
 
-  //get nFaces
-  while(temp.find("nFaces") == vtkstd::string::npos)
-    {
-    tempStringStruct = this->GetLine(input);
-    temp = tempStringStruct->value;
-    delete tempStringStruct;
-    }
-  temp.erase(temp.begin()+temp.find(";")); //remove ;
-  tokenizer << temp;
-  //while(tokenizer >> token);
-  while(!tokenizer.eof())
-    {
-    tokenizer >> token;
-    }
-  tokenizer.clear();
-  tokenizer << token;
-  tokenizer >> nFaces;
+    // Create a boundary mesh only if selected in the list for display
+    vtkFoamEntry& boundaryEntryI = boundaryDict.entry(boundaryIndex);
+    const vtkStdString& boundaryNameI = boundaryEntryI.keyword();
+    if(!this->GetPatchArrayStatus(boundaryNameI.c_str()))
+      {
+      continue;
+      }
 
-  //get startface
-  tempStringStruct = this->GetLine(input);
-  temp = tempStringStruct->value;
-  delete tempStringStruct;
+    int startFace;
+    startFaceEntry >> startFace;
 
-  //look for "startFaces"
-  while(temp.find("startFace") == vtkstd::string::npos)
-    {
-    tempStringStruct = this->GetLine(input);
-    temp = tempStringStruct->value;
-    delete tempStringStruct;
-    }
-  temp.erase(temp.begin()+temp.find(";")); //remove ;
-  tokenizer.clear();
-  tokenizer << temp;
-  //while(tokenizer >> token);
-  while(!tokenizer.eof())
-    {
-    tokenizer >> token;
-    }
-  tokenizer.clear();
-  tokenizer << token;
-  tokenizer >> this->StartFace;
+    int endFace = startFace + nFaces;
 
-  //Create the mesh
-  int j, k;
-  vtkTriangle * triangle;
-  vtkQuad * quad;
-  vtkPolygon * polygon;
-  int endFace = this->StartFace + nFaces;
-  //loop through each face
-  for(j = this->StartFace; j < endFace; j++)
-    {
+    //Create the mesh
+    boundaryMesh->value.push_back(vtkUnstructuredGrid::New());
+    vtkUnstructuredGrid& bm = *boundaryMesh->value.back();
+    bm.Allocate(nFaces);
 
-    //triangle
-    if(this->FacePoints->value[j].size() == 3)
+    // aloocate array for converting int vector to vtkIdType vector:
+    // workaround for 64bit machines
+    int maxNFacePoints = 0;
+    for(int j = startFace; j < endFace; j++)
       {
-      triangle = vtkTriangle::New();
-      for(k = 0; k < 3; k++)
+      const int nFacePoints = facesPoints->value[j].size();
+      if(nFacePoints > maxNFacePoints)
         {
-        triangle->GetPointIds()->SetId(k, this->FacePoints->value[j][k]);
+        maxNFacePoints = nFacePoints;
         }
-      boundaryMesh->InsertNextCell(triangle->GetCellType(),
-      triangle->GetPointIds());
-      triangle->Delete();
       }
+    vtkIdType* facePointsVtkId = new vtkIdType[maxNFacePoints];
 
-    //quad
-    else if(this->FacePoints->value[j].size() == 4)
+    //loop through each face
+    for(int j = startFace; j < endFace; j++)
       {
-      quad = vtkQuad::New();
-      for(k = 0; k < 4; k++)
+      vtkstd::vector<int>& facePoints = facesPoints->value[j];
+      size_t nFacePoints = facePoints.size();
+
+      for(size_t k = 0; k < nFacePoints; k++)
         {
-        quad->GetPointIds()->SetId(k, this->FacePoints->value[j][k]);
+        facePointsVtkId[k] = facePoints[k];
         }
-      boundaryMesh->InsertNextCell(quad->GetCellType(),
-      quad->GetPointIds());
-      quad->Delete();
-      }
 
-    //polygon
-    else
-      {
-      polygon = vtkPolygon::New();
-      for(k = 0; k < (int)this->FacePoints->value[j].size(); k++)
+      //triangle
+      if(nFacePoints == 3)
+        {
+        bm.InsertNextCell(VTK_TRIANGLE, 3, facePointsVtkId);
+        }
+      // quad
+      else if(nFacePoints == 4)
         {
-        polygon->GetPointIds()->InsertId(k, this->FacePoints->value[j][k]);
+        bm.InsertNextCell(VTK_QUAD, 4, facePointsVtkId);
+        }
+      //polygon
+      else
+        {
+        bm.InsertNextCell(VTK_POLYGON, nFacePoints,
+          facePointsVtkId);
         }
-      boundaryMesh->InsertNextCell(polygon->GetCellType(),
-      polygon->GetPointIds());
-      polygon->Delete();
       }
-    }
+    delete [] facePointsVtkId;
 
-  //set points for boundary
-  boundaryMesh->SetPoints(this->Points);
-  //add size of mesh
-  this->SizeOfBoundary->value.push_back(boundaryMesh->GetNumberOfCells());
-  input->close();
-  delete input;
+    //set points for boundary
+    bm.SetPoints(points);
+    }
   vtkDebugMacro(<<"Boundary mesh created");
   return boundaryMesh;
 }
@@ -2521,148 +4460,69 @@
 //  returns a requested point zone mesh
 //
 // ****************************************************************************
-vtkUnstructuredGrid * vtkOpenFOAMReader::GetPointZoneMesh(int timeState,
-                                                          int pointZoneIndex)
+bool vtkOpenFOAMReader::GetPointZoneMesh(unstructuredGridVector* pointZoneMesh,
+  vtkPoints* points, int timeState)
 {
-  vtkUnstructuredGrid * pointZoneMesh = vtkUnstructuredGrid::New();
-  vtkstd::string pointZonesPath = this->PathPrefix->value[timeState]+
-                               "/polyMesh/pointZones";
-  vtkDebugMacro(<<"Create point zone mesh: "<<pointZonesPath.c_str());
+  vtkDebugMacro(<<"Create point zone mesh");
 
-  vtkstd::string temp;
-  stdString* tempStringStruct;
-  bool binaryWriteFormat;
-  ifstream * input = new ifstream(pointZonesPath.c_str(), ios::in VTK_IOS_NOCREATE);
-  //make sure file exists
-  if(input->fail())
-    {
-    input->close();
-    delete input;
-    return pointZoneMesh;
-    }
+  vtkFoamDict* pointZoneDictPtr
+    = this->GatherBlocks("pointZones", timeState, false);
 
-  //determine if file is binary or ascii
-  while(temp.find("format") == vtkstd::string::npos)
+  if(pointZoneDictPtr == NULL)
     {
-    tempStringStruct = this->GetLine(input);
-    temp = tempStringStruct->value;
-    delete tempStringStruct;
+    // not an error
+    return true;
     }
-  input->close();
 
-  //reopen file in correct format
-  if(temp.find("binary") != vtkstd::string::npos)
-    {
-#ifdef _WIN32
-    input->open(pointZonesPath.c_str(), ios::binary | ios::in VTK_IOS_NOCREATE);
-#else
-    input->open(pointZonesPath.c_str(), ios::in VTK_IOS_NOCREATE);
-#endif
-    binaryWriteFormat = true;
-    }
-  else
-    {
-    input->open(pointZonesPath.c_str(),ios::in);
-    binaryWriteFormat = false;
-    }
+  vtkFoamDict& pointZoneDict = *pointZoneDictPtr;
+  size_t nPointZones = pointZoneDict.size();
 
-  vtkstd::string token;
-  vtksys_ios::stringstream tokenizer;
-  vtkVertex * pointCell;
-  int tempElement;
-  vtkstd::vector< vtkstd::vector < int > > tempElementZones;
-  int numElement;
-
-  //find desired mesh entry
-  while(temp.find(this->PointZoneNames->value[pointZoneIndex]) == 
-        vtkstd::string::npos)
-    {
-    tempStringStruct = this->GetLine(input);
-    temp = tempStringStruct->value;
-    delete tempStringStruct;
-    }
-  tempStringStruct = this->GetLine(input);
-  temp = tempStringStruct->value;
-  delete tempStringStruct;//throw out {
-
-  tempStringStruct = this->GetLine(input);
-  temp = tempStringStruct->value;
-  delete tempStringStruct;//type
-
-  tempStringStruct = this->GetLine(input);
-  temp = tempStringStruct->value;
-  delete tempStringStruct;//label
-
-  tempStringStruct = this->GetLine(input);
-  temp = tempStringStruct->value;
-  delete tempStringStruct;//number of elements or {
-
-  //number of elements
-  if(temp.find("}") == vtkstd::string::npos)
-    {
-    tokenizer << temp;
-    tokenizer >> numElement;
-    if(numElement == 0)
+  for(size_t i = 0; i < nPointZones; i++)
+    {
+    // look up point labels
+    vtkFoamDict& dict = pointZoneDict.entry(i).dictionary();
+    vtkFoamEntry& pointLabelsEntry = dict.lookup("pointLabels");
+    if(!pointLabelsEntry.found())
       {
-      input->close();
-      delete input;
-      return NULL;
+      vtkErrorMacro(<<"pointLabels not found in pointZones");
+      return false;
       }
-
-    //binary data
-    if(binaryWriteFormat)
+    // skip if the list is empty
+    if(pointLabelsEntry.firstValue().type() == vtkFoamEntryValue::EMPTYLIST)
       {
-      input->get(); //parenthesis
-      for(int j = 0; j < numElement; j++)
-        {
-        input->read((char *) &tempElement, sizeof(int));
-        pointCell = vtkVertex::New();
-        pointCell->GetPointIds()->SetId(0,tempElement);
-        pointZoneMesh->InsertNextCell(pointCell->GetCellType(),
-                                      pointCell->GetPointIds());
-        pointCell->Delete();
-        }
+      continue;
       }
-
-    //ascii data
-    else
+    if(pointLabelsEntry.firstValue().type() != vtkFoamEntryValue::LABELLIST)
       {
-      tempStringStruct = this->GetLine(input);
-      temp = tempStringStruct->value;
-      delete tempStringStruct;//THROW OUT (
+      vtkErrorMacro(<<"pointLabels not of type labelList: type = "
+        << pointLabelsEntry.firstValue().type());
+      return false;
+      }
 
-      //GET EACH ELEMENT & ADD TO VECTOR
-      for(int j = 0; j < numElement; j++)
-        {
-        tempStringStruct = this->GetLine(input);
-        temp = tempStringStruct->value;
-        delete tempStringStruct;
+    const vtkstd::vector<int>& labels = pointLabelsEntry.labelList();
 
-        tokenizer.clear();
-        tokenizer << temp;
-        tokenizer >> tempElement;
-        pointCell = vtkVertex::New();
-        pointCell->GetPointIds()->SetId(0,tempElement);
-        pointZoneMesh->InsertNextCell(pointCell->GetCellType(),
-                                      pointCell->GetPointIds());
-        pointCell->Delete();
-        }
-      }
-    }
+    // allocate new grid: we do not use resize() beforehand since it
+    // could lead to undefined pointer if we return by error
+    pointZoneMesh->value.push_back(vtkUnstructuredGrid::New());
+    vtkUnstructuredGrid& pzm = *pointZoneMesh->value.back();
 
-  //there is no entry
-  else
-    {
-    input->close();
-    delete input;
-    return NULL;
+    // set pointZone size
+    size_t nPoints = labels.size();
+    pzm.Allocate(nPoints);
+
+    // insert points
+    for(size_t j = 0; j < nPoints; j++)
+      {
+      vtkIdType pointLabel = labels[j];
+      pzm.InsertNextCell(VTK_VERTEX, 1, &pointLabel);
+      }
+    pzm.SetPoints(points);
     }
-  //set point zone points
-  pointZoneMesh->SetPoints(Points);
-  input->close();
-  delete input;
+
+  delete pointZoneDictPtr;
+
   vtkDebugMacro(<<"Point zone mesh created");
-  return pointZoneMesh;
+  return true;
 }
 
 // ****************************************************************************
@@ -2672,189 +4532,108 @@
 //  returns a requested face zone mesh
 //
 // ****************************************************************************
-vtkUnstructuredGrid * vtkOpenFOAMReader::GetFaceZoneMesh(int timeState,
-                                                         int faceZoneIndex)
+bool vtkOpenFOAMReader::GetFaceZoneMesh(
+  unstructuredGridVector* faceZoneMesh, intVectorVector* facesPoints,
+  vtkPoints* points, int timeState)
 {
-  vtkUnstructuredGrid * faceZoneMesh = vtkUnstructuredGrid::New();
-  vtkstd::string faceZonesPath = this->PathPrefix->value +
-                                 this->PolyMeshFacesDir->value[timeState] +
-                                 "/polyMesh/faceZones";
-  vtkDebugMacro(<<"Create face zone mesh: "<<faceZonesPath.c_str());
-
-  vtkstd::string temp;
-  stdString* tempStringStruct;
-  bool binaryWriteFormat;
-  ifstream * input = new ifstream(faceZonesPath.c_str(), ios::in VTK_IOS_NOCREATE);
-  //make sure file exists
-  if(input->fail())
-    {
-    input->close();
-    delete input;
-    return faceZoneMesh;
-    }
-
-  //determine if file is binary or ascii
-  while(temp.find("format") == vtkstd::string::npos)
-    {
-    tempStringStruct = this->GetLine(input);
-    temp = tempStringStruct->value;
-    delete tempStringStruct;
-    }
-  input->close();
-
-  //reopen file in correct format
-  if(temp.find("binary") != vtkstd::string::npos)
-    {
-#ifdef _WIN32
-    input->open(faceZonesPath.c_str(), ios::binary | ios::in VTK_IOS_NOCREATE);
-#else
-    input->open(faceZonesPath.c_str(), ios::in);
-#endif
-    binaryWriteFormat = true;
-    }
-  else
-    {
-    input->open(faceZonesPath.c_str(),ios::in);
-    binaryWriteFormat = false;
-    }
+  vtkDebugMacro(<<"Create face zone mesh");
 
-  vtkstd::string token;
-  vtksys_ios::stringstream tokenizer;
-  vtkstd::vector< int > faceZone;
-  int tempElement;
-  vtkstd::vector< vtkstd::vector < int > > tempElementZones;
-  int numElement;
+  vtkFoamDict* faceZoneDictPtr
+    = this->GatherBlocks("faceZones", timeState, false);
 
-  //find desired mesh entry
-  while(temp.find(this->FaceZoneNames->value[faceZoneIndex]) == 
-        vtkstd::string::npos)
+  if(faceZoneDictPtr == NULL)
     {
-    tempStringStruct = this->GetLine(input);
-    temp = tempStringStruct->value;
-    delete tempStringStruct;
+    // not an error
+    return true;
     }
 
-  tempStringStruct = this->GetLine(input);
-  temp = tempStringStruct->value;
-  delete tempStringStruct;//throw out {
+  vtkFoamDict& faceZoneDict = *faceZoneDictPtr;
+  size_t nFaceZones = faceZoneDict.size();
 
-  tempStringStruct = this->GetLine(input);
-  temp = tempStringStruct->value;
-  delete tempStringStruct;//type
-
-  tempStringStruct = this->GetLine(input);
-  temp = tempStringStruct->value;
-  delete tempStringStruct;//label
-
-  tempStringStruct = this->GetLine(input);
-  temp = tempStringStruct->value;
-  delete tempStringStruct;//number of values or flipmap
-
-  if(temp.find("flipMap") == vtkstd::string::npos)
+  for(size_t i = 0; i < nFaceZones; i++)
     {
-    //number of elements
-    tokenizer << temp;
-    tokenizer >> numElement;
-    if(numElement == 0)
+    // look up face labels
+    vtkFoamDict& dict = faceZoneDict.entry(i).dictionary();
+    vtkFoamEntry& faceLabelsEntry = dict.lookup("faceLabels");
+    if(!faceLabelsEntry.found())
       {
-      input->close();
-      delete input;
-      return NULL;
+      delete faceZoneDictPtr;
+      vtkErrorMacro(<<"faceLabels not found in faceZones");
+      return false;
       }
-
-    //binary
-    if(binaryWriteFormat)
+    // skip if the list is empty
+    if(faceLabelsEntry.firstValue().type() == vtkFoamEntryValue::EMPTYLIST)
       {
-      input->get(); //parenthesis
-      for(int j = 0; j < numElement; j++)
-        {
-        input->read((char *) &tempElement, sizeof(int));
-        faceZone.push_back(tempElement);
-        }
+      continue;
       }
-
-    //ascii
-    else
+    if(faceLabelsEntry.firstValue().type() != vtkFoamEntryValue::LABELLIST)
       {
-      //THROW OUT (
-      tempStringStruct = this->GetLine(input);
-      temp = tempStringStruct->value;
-      delete tempStringStruct;
-
-      //get each element & add to vector
-      for(int j = 0; j < numElement; j++)
-        {
-        tempStringStruct = this->GetLine(input);
-        temp = tempStringStruct->value;
-        delete tempStringStruct;
-
-        tokenizer.clear();
-        tokenizer << temp;
-        tokenizer >> tempElement;
-        faceZone.push_back(tempElement);
-        }
+      delete faceZoneDictPtr;
+      vtkErrorMacro(<<"faceLabels not of type labelList");
+      return false;
       }
-    }
 
-  //Create the mesh
-  int k;
-  vtkTriangle * triangle;
-  vtkQuad * quad;
-  vtkPolygon * polygon;
+    const vtkstd::vector<int>& labels = faceLabelsEntry.labelList();
 
-  //LOOP THROUGH EACH FACE
-  for(int j = 0; j < (int)faceZone.size(); j++)
-    {
+    // allocate new grid: we do not use resize() beforehand since it
+    // could lead to undefined pointer if we return by error
+    faceZoneMesh->value.push_back(vtkUnstructuredGrid::New());
+    vtkUnstructuredGrid& fzm = *faceZoneMesh->value.back();
+
+    // set faceZone size
+    size_t nFaces = labels.size();
+    fzm.Allocate(nFaces);
 
-    //Triangular Face
-    if(this->FacePoints->value[faceZone[j]].size() == 3)
+    // aloocate array for converting int vector to vtkIdType vector:
+    // workaround for 64bit machines
+    int maxNFacePoints = 0;
+    for(size_t j = 0; j < nFaces; j++)
       {
-      triangle = vtkTriangle::New();
-      for(k = 0; k < 3; k++)
+      const int nFacePoints = facesPoints->value[labels[j]].size();
+      if(nFacePoints > maxNFacePoints)
         {
-        triangle->GetPointIds()->SetId(k, this->FacePoints->value[
-        faceZone[j]][k]);
+        maxNFacePoints = nFacePoints;
         }
-      faceZoneMesh->InsertNextCell(triangle->GetCellType(),
-        triangle->GetPointIds());
-      triangle->Delete();
       }
+    vtkIdType* facePointsVtkId = new vtkIdType[maxNFacePoints];
 
-    //Quadraic Face
-    else if(this->FacePoints->value[faceZone[j]].size() == 4)
+    // insert faces
+    for(size_t j = 0; j < nFaces; j++)
       {
-      quad = vtkQuad::New();
-      for(k = 0; k < 4; k++)
+      vtkstd::vector<int>& facePoints = facesPoints->value[labels[j]];
+      size_t nFacePoints = facePoints.size();
+
+      for(size_t k = 0; k < nFacePoints; k++)
         {
-        quad->GetPointIds()->SetId(k,
-          this->FacePoints->value[faceZone[j]][k]);
+        facePointsVtkId[k] = facePoints[k];
         }
-      faceZoneMesh->InsertNextCell(quad->GetCellType(),
-        quad->GetPointIds());
-      quad->Delete();
-      }
 
-    //Polygonal Face
-    else
-      {
-      polygon = vtkPolygon::New();
-      for(k = 0; k < (int)this->FacePoints->value[faceZone[j]].size(); k++)
+      //triangle
+      if(nFacePoints == 3)
+        {
+        fzm.InsertNextCell(VTK_TRIANGLE, 3, facePointsVtkId);
+        }
+      // quad
+      else if(nFacePoints == 4)
+        {
+        fzm.InsertNextCell(VTK_QUAD, 4, facePointsVtkId);
+        }
+      //polygon
+      else
         {
-        polygon->GetPointIds()->InsertId(k, this->FacePoints->value[
-          faceZone[j]][k]);
+        fzm.InsertNextCell(VTK_POLYGON, nFacePoints,
+          facePointsVtkId);
         }
-      faceZoneMesh->InsertNextCell(polygon->GetCellType(),
-        polygon->GetPointIds());
-      polygon->Delete();
       }
+
+    delete [] facePointsVtkId;
+    fzm.SetPoints(points);
     }
 
-  //set the face zone points
-  faceZoneMesh->SetPoints(this->Points);
-  input->close();
-  delete input;
+  delete faceZoneDictPtr;
+
   vtkDebugMacro(<<"Face zone mesh created");
-  return faceZoneMesh;
+  return true;
 }
 
 // ****************************************************************************
@@ -2864,628 +4643,590 @@
 //  returns a requested cell zone mesh
 //
 // ****************************************************************************
-vtkUnstructuredGrid * vtkOpenFOAMReader::GetCellZoneMesh(int timeState,
-                                                         int cellZoneIndex)
+bool vtkOpenFOAMReader::GetCellZoneMesh(
+  unstructuredGridVector* cellZoneMesh, faceVectorVector* facesOfCell,
+  intVectorVector* facesPoints, vtkPoints* points, int timeState)
 {
-  vtkUnstructuredGrid * cellZoneMesh = vtkUnstructuredGrid::New();
-  vtkstd::string cellZonesPath = this->PathPrefix->value +
-                                 this->PolyMeshFacesDir->value[timeState] +
-                                 "/polyMesh/cellZones";
-  vtkDebugMacro(<<"Create cell zone mesh: "<<cellZonesPath.c_str());
-  vtkstd::string temp;
-  stdString* tempStringStruct;
-  bool binaryWriteFormat;
-  ifstream * input = new ifstream(cellZonesPath.c_str(), ios::in VTK_IOS_NOCREATE);
-  //make sure file exists
-  if(input->fail())
-    {
-    input->close();
-    delete input;
-    return cellZoneMesh;
-    }
-
-  //determine if file is binary or ascii
-  while(temp.find("format") == vtkstd::string::npos)
-    {
-    tempStringStruct = this->GetLine(input);
-    temp = tempStringStruct->value;
-    delete tempStringStruct;
-    }
-  input->close();
+  vtkDebugMacro(<<"Create cell zone mesh");
 
-  //reopen file in correct format
-  if(temp.find("binary") != vtkstd::string::npos)
+  vtkFoamDict* cellZoneDictPtr
+    = this->GatherBlocks("cellZones", timeState, false);
+
+  if(cellZoneDictPtr == NULL)
     {
-#ifdef _WIN32
-    input->open(cellZonesPath.c_str(), ios::binary | ios::in VTK_IOS_NOCREATE);
-#else
-    input->open(cellZonesPath.c_str(), ios::in VTK_IOS_NOCREATE);
-#endif
-    binaryWriteFormat = true;
+    // not an error
+    return true;
     }
-  else
+
+  vtkFoamDict& cellZoneDict = *cellZoneDictPtr;
+  size_t nCellZones = cellZoneDict.size();
+
+  for(size_t i = 0; i < nCellZones; i++)
     {
-    input->open(cellZonesPath.c_str(),ios::in);
-    binaryWriteFormat = false;
+    // look up cell labels
+    vtkFoamDict& dict = cellZoneDict.entry(i).dictionary();
+    vtkFoamEntry& cellLabelsEntry = dict.lookup("cellLabels");
+    if(!cellLabelsEntry.found())
+      {
+      delete cellZoneDictPtr;
+      vtkErrorMacro(<<"cellLabels not found in cellZones");
+      return false;
+      }
+    // skip if the list is empty
+    if(cellLabelsEntry.firstValue().type() == vtkFoamEntryValue::EMPTYLIST)
+      {
+      continue;
+      }
+    if(cellLabelsEntry.firstValue().type() != vtkFoamEntryValue::LABELLIST)
+      {
+      delete cellZoneDictPtr;
+      vtkErrorMacro(<<"cellLabels not of type labelList");
+      return false;
+      }
+
+    const vtkstd::vector<int>& labels = cellLabelsEntry.labelList();
+
+    // allocate new grid: we do not use resize() beforehand since it
+    // could lead to undefined pointer if we return by error
+    cellZoneMesh->value.push_back(vtkUnstructuredGrid::New());
+    vtkUnstructuredGrid* czm = cellZoneMesh->value.back();
+
+    // set cellZone size
+    size_t nCells = labels.size();
+    czm->Allocate(nCells);
+
+    // insert cells
+    for(size_t j = 0; j < nCells; j++)
+      {
+      this->InsertCellToGrid(czm, labels[j], facesOfCell, facesPoints);
+      }
+
+    //set cell zone points
+    cellZoneMesh->value[i]->SetPoints(points);
     }
 
-  vtkstd::string token;
-  vtksys_ios::stringstream tokenizer;
-  vtkstd::vector< int > cellZone;
-  int tempElement;
-  vtkstd::vector< vtkstd::vector < int > > tempElementZones;
-  int numElement;
+  delete cellZoneDictPtr;
+  vtkDebugMacro(<<"Cell zone mesh created");
+  return true;
+}
+
+vtkPolyData* vtkOpenFOAMReader::MakeLagrangianMesh(int timeState)
+{
+  vtkStdString positionsPath = *this->PathPrefix
+    + this->TimeNames->GetValue(timeState) +  "/lagrangian/positions";
 
-  //find desired mesh entry
-  while(temp.find(this->CellZoneNames->value[cellZoneIndex]) == 
-        vtkstd::string::npos)
+  vtkFoamIOobject io;
+  if(!(io.open(positionsPath) || io.open(positionsPath + ".gz")))
+    {
+    vtkErrorMacro(<<"Error opening " << io.fileName().c_str() << ": "
+      << io.error().str().c_str());
+    return NULL;
+    }
+  vtkFoamDict dict;
+  if(!dict.read(io))
+    {
+    vtkErrorMacro(<<"Error reading line " << io.lineNumber()
+      << " of " << io.fileName().c_str() << ": " << io.error().str().c_str());
+    return NULL;
+    }
+  if(dict.type() != vtkFoamDict::VECTORLIST)
     {
-    tempStringStruct = this->GetLine(input);
-    temp = tempStringStruct->value;
-    delete tempStringStruct;
+    vtkErrorMacro(<<"The file type of " << io.fileName().c_str()
+      << " is not a vectorList");
+    return NULL;
     }
-  tempStringStruct = this->GetLine(input);
-  temp = tempStringStruct->value;
-  delete tempStringStruct;//throw out {
 
-  tempStringStruct = this->GetLine(input);
-  temp = tempStringStruct->value;
-  delete tempStringStruct;//type
+  vtkDoubleArray& vl = dict.vectorList();
+  const vtkIdType nParticles = vl.GetNumberOfTuples();
+  //instantiate the points class
+  vtkPoints* points = vtkPoints::New();
+  points->SetNumberOfPoints(nParticles);
+  vtkPolyData* lagrangianMesh = vtkPolyData::New();
+  lagrangianMesh->Allocate(nParticles);
 
-  tempStringStruct = this->GetLine(input);
-  temp = tempStringStruct->value;
-  delete tempStringStruct;
+  for(vtkIdType i = 0, index = 0; i < nParticles; i++, index += 3)
+    {
+    points->SetPoint(i, vl.GetPointer(index));
+    lagrangianMesh->InsertNextCell(VTK_VERTEX, 1, &i);
+    }
 
-  tempStringStruct = this->GetLine(input);
-  temp = tempStringStruct->value;
-  delete tempStringStruct;
+  lagrangianMesh->SetPoints(points);
+  points->Delete();
+  return lagrangianMesh;
+}
 
-  //number of elements
-  tokenizer << temp;
-  tokenizer >> numElement;
+bool vtkOpenFOAMReader::GetLagrangianVariableAtTimestep(
+  vtkPolyData* lagrangianMesh, const char* varNameIn, int timeState)
+{
+  vtkStdString varPath = *this->PathPrefix
+    + this->TimeNames->GetValue(timeState) + "/lagrangian/"
+    + vtkStdString(varNameIn);
 
-  //binary
-  if(binaryWriteFormat)
+  // open the file
+  vtkFoamIOobject io;
+  if(!io.open(varPath))
     {
-    input->get(); //parenthesis
-    for(int j = 0; j < numElement; j++)
-      {
-      input->read((char *) &tempElement, sizeof(int));
-      cellZone.push_back(tempElement);
-      }
+    vtkErrorMacro(<<"Error opening " << io.fileName().c_str() << ": "
+      << io.error().str().c_str());
+    return false;
     }
 
-  //ascii
-  else
+  // if the variable is disabled on selection panel then skip it
+  vtkStdString selectionName
+    = vtkStdString("lagrangian/") + io.objectName();
+  if(this->PointDataArraySelection->ArrayExists(selectionName.c_str())
+    && !this->PointDataArraySelection->ArrayIsEnabled(selectionName.c_str()))
     {
-    tempStringStruct = this->GetLine(input);
-    temp = tempStringStruct->value;
-    delete tempStringStruct;//throw out (
-
-    //get each element & add to vector
-    for(int j = 0; j < numElement; j++)
-      {
-      tempStringStruct = this->GetLine(input);
-      temp = tempStringStruct->value;
-      delete tempStringStruct;
-
-      tokenizer.clear();
-      tokenizer << temp;
-      tokenizer >> tempElement;
-      cellZone.push_back(tempElement);
-      }
+    return true;
     }
 
-  //Create the mesh
-  bool foundDup = false;
-  vtkstd::vector< int > cellPoints;
-  vtkstd::vector< int > tempFaces[2];
-  vtkstd::vector< int > firstFace;
-  int pivotPoint = 0;
-  int i, j, k, l, pCount;
-  int faceCount = 0;
+  // read the field file into dictionary
+  vtkFoamDict dict;
+  if(!dict.read(io))
+    {
+    vtkErrorMacro(<<"Error reading line " << io.lineNumber()
+      << " of " << io.fileName().c_str() << ": " << io.error().str().c_str());
+    return false;
+    }
 
-  //Create Mesh
-  for(i = 0; i < (int)cellZone.size(); i++)  //each cell
+  // set lagrangian values
+  if(dict.type() != vtkFoamDict::SCALARLIST
+    && dict.type() != vtkFoamDict::VECTORLIST)
     {
-    //calculate total points for all faces of a cell
-    //used to determine cell type
-    int totalPointCount = 0;
-    for(j = 0; j < (int)this->FacesOfCell->value[
-      cellZone[i]].size(); j++)  //each face
-      {
-      totalPointCount += (int)this->FacePoints->value[this->FacesOfCell->value[
-        cellZone[i]][j].faceIndex].size();
-      }
+    vtkErrorMacro(<<"Unsupported lagrangian field type");
+    return false;
+    }
 
-      // using cell type - order points, create cell, add to mesh
+  vtkDoubleArray* lData = (vtkDoubleArray *)dict.ptr();
 
-    //OFhex | vtkHexahedron
-    if (totalPointCount == 24)
-      {
-      faceCount = 0;
+  // GetNumberOfTuples() works for both scalar and vector
+  const int nParticles = lData->GetNumberOfTuples();
+  if(nParticles != lagrangianMesh->GetNumberOfCells())
+    {
+    vtkErrorMacro(<<"Sizes of lagrangian mesh and field don't match");
+    return false;
+    }
 
-      //get first face
-      for(j = 0; j < (int)this->FacePoints->value[this->FacesOfCell->value[
-        cellZone[i]][0].faceIndex].size(); j++)
-        {
-        firstFace.push_back(this->FacePoints->value[this->FacesOfCell->value[
-          cellZone[i]][0].faceIndex][j]);
-        }
+  lData->SetName(selectionName.c_str());
+  lagrangianMesh->GetPointData()->AddArray(lData);
+  lData->Delete();
+
+  return true;
+}
 
-      //-if it is a neighbor face flip it
-      if(FacesOfCell->value[i][0].neighborFace)
+// return 0 if there's any error, 1 if success
+int vtkOpenFOAMReader::CreateDataSet(vtkMultiBlockDataSet *output,
+  int timeState)
+{
+  int oldTimeState = this->OldTimeStep;
+  vtkPolyData* lagrangianMesh = NULL;
+  bool updateVariables = true;
+
+  // determine if we need to reconstruct meshes
+  if(this->CacheMesh && oldTimeState >= 0
+    && this->PolyMeshFacesDir->value[timeState]
+       == this->PolyMeshFacesDir->value[oldTimeState]
+    && this->FaceOwner != NULL
+    && this->PatchSelectionStatus == this->PatchSelectionOldStatus)
+    {
+    if((timeState == oldTimeState) 
+       && (this->CellSelectionStatus == this->CellSelectionOldStatus)
+       && (this->PointSelectionStatus == this->PointSelectionOldStatus))
+      {
+      updateVariables = false;
+      }
+    else
+      {
+      // clean up arrays of the previous timestep
+      // Check if Internal Mesh Exists first...
+      if(this->InternalMesh != NULL)
         {
-        int tempPop;
-        for(k = 0; k < (int)firstFace.size() - 1; k++)
+        vtkCellData* internalCellData = this->InternalMesh->GetCellData();
+        while(internalCellData->GetArrayName(0))
           {
-          tempPop = firstFace[firstFace.size()-1];
-          firstFace.pop_back();
-          firstFace.insert(firstFace.begin()+1+k, tempPop);
+          internalCellData->RemoveArray(internalCellData->GetArrayName(0));
           }
         }
-
-      //add first face to cell points
-      for(j =0; j < (int)firstFace.size(); j++)
-        {
-        cellPoints.push_back(firstFace[j]);
-        }
-
-      for(int pointCount = 0;pointCount < (int)firstFace.size();pointCount++)
+      // Check if Boundary Mesh Exists first...
+      if(this->BoundaryMesh != NULL)
         {
-        //find the 2 other faces containing each point - start with face 1
-        for(j = 1; j < (int)this->FacesOfCell->value[
-          cellZone[i]].size(); j++)  //each face
+        for(size_t i = 0; i < this->BoundaryMesh->value.size(); i++)
           {
-          for(k = 0; k < (int)this->FacePoints->value[this->FacesOfCell->
-              value[cellZone[i]][j].faceIndex].size(); k++)
+          vtkCellData* boundaryCellData
+            = this->BoundaryMesh->value[i]->GetCellData();
+          while(boundaryCellData->GetArrayName(0))
             {
-            if(firstFace[pointCount] == this->FacePoints->
-               value[this->FacesOfCell->value[cellZone[i]][j].faceIndex][k])
-              {
-              //another face with the point
-              for(l = 0; l < (int)this->FacePoints->value[this->FacesOfCell->
-                                  value[cellZone[i]][j].faceIndex].size(); l++)
-                {
-                tempFaces[faceCount].push_back(this->FacePoints->value[
-                  this->FacesOfCell->value[cellZone[i]][j].faceIndex]
-                  [l]);
-                }
-              faceCount++;
-              }
+            boundaryCellData->RemoveArray(boundaryCellData->GetArrayName(0));
             }
           }
-
-        //locate the pivot point contained in faces 0 & 1
-        for(j = 0; j < (int)tempFaces[0].size(); j++)
+        }
+      // if only point coordinates change refresh point vector
+      if(this->PolyMeshPointsDir->value[timeState]
+        != this->PolyMeshPointsDir->value[oldTimeState])
+        {
+        //get the points
+        vtkPoints* points = this->GetPoints(timeState);
+        if(points == NULL)
           {
-          for(k = 0; k < (int)tempFaces[1].size(); k++)
-            {
-            if(tempFaces[0][j] == tempFaces[1][k] && tempFaces[0][j] !=
-              firstFace[pointCount])
-              {
-              pivotPoint = tempFaces[0][j];
-              break;
-              }
-            }
+          return 0;
           }
-        cellPoints.push_back(pivotPoint);
-        tempFaces[0].clear();
-        tempFaces[1].clear();
-        faceCount=0;
-        }
 
-        //create the hex cell and insert it into the mesh
-        vtkHexahedron * hexahedron = vtkHexahedron::New();
-        for(pCount = 0; pCount < (int)cellPoints.size(); pCount++)
+        // refresh the points in each mesh
+        // Check if Internal Mesh exists first....
+        if(this->InternalMesh != NULL)
           {
-          hexahedron->GetPointIds()->SetId(pCount, cellPoints[pCount]);
+          this->InternalMesh->SetPoints(points);
           }
-        cellZoneMesh->InsertNextCell(hexahedron->GetCellType(),
-                                     hexahedron->GetPointIds());
-        hexahedron->Delete();
-        cellPoints.clear();
-        firstFace.clear();
-        }
-
-    //OFprism | vtkWedge
-    else if (totalPointCount == 18)
-      {
-      faceCount = 0;
-      int index = 0;
-
-      //find first triangular face
-      for(j = 0; j < (int)this->FacesOfCell->value[
-        cellZone[i]].size(); j++)  //each face
-        {
-        if((int)this->FacePoints->value[this->FacesOfCell->value[
-          cellZone[i]][j].faceIndex].size() == 3)
+        // Check if Boundary Mesh exists first....
+        if(this->BoundaryMesh != NULL)
           {
-          for(k = 0; k < (int)this->FacePoints->value[this->FacesOfCell->value[
-            cellZone[i]][j].faceIndex].size(); k++)
+          for(size_t i = 0; i < this->BoundaryMesh->value.size(); i++)
             {
-            firstFace.push_back(this->FacePoints->value[this->FacesOfCell->
-                                  value[cellZone[i]][j].faceIndex][k]);
-            index = j;
+            this->BoundaryMesh->value[i]->SetPoints(points);
             }
-          break;
           }
-        }
-
-      //-if it is a neighbor face flip it
-      if(this->FacesOfCell->value[i][0].neighborFace)
-        {
-        int tempPop;
-        for(k = 0; k < (int)firstFace.size() - 1; k++)
+        if(this->PointZoneMesh != NULL)
           {
-          tempPop = firstFace[firstFace.size()-1];
-          firstFace.pop_back();
-          firstFace.insert(firstFace.begin()+1+k, tempPop);
+          for(size_t i = 0; i < this->PointZoneMesh->value.size(); i++)
+            {
+            this->PointZoneMesh->value[i]->SetPoints(points);
+            }
           }
-        }
-
-      //add first face to cell points
-      for(j =0; j < (int)firstFace.size(); j++)
-        {
-        cellPoints.push_back(firstFace[j]);
-        }
-
-      for(int pointCount = 0;pointCount < (int)firstFace.size();pointCount++)
-        {
-        //find the 2 other faces containing each point
-        for(j = 0; j < (int)this->FacesOfCell->value[
-          cellZone[i]].size(); j++)  //each face
+        if(this->FaceZoneMesh != NULL)
           {
-          for(k = 0; k < (int)this->FacePoints->value[this->FacesOfCell->value[
-            cellZone[i]][j].faceIndex].size(); k++) //each point
+          for(size_t i = 0; i < this->FaceZoneMesh->value.size(); i++)
             {
-            if(firstFace[pointCount] == this->FacePoints->
-               value[this->FacesOfCell->value[cellZone[i]][j].faceIndex][k] &&
-               j != index)
-              {
-              //another face with point
-              for(l = 0; l < (int)this->FacePoints->value[this->
-                  FacesOfCell->value[cellZone[i]][j].faceIndex].size(); l++)
-                {
-                tempFaces[faceCount].push_back(this->FacePoints->value[
-                  this->FacesOfCell->value[cellZone[i]][j].faceIndex]
-                  [l]);
-                }
-              faceCount++;
-              }
+            this->FaceZoneMesh->value[i]->SetPoints(points);
             }
           }
-
-        //locate the pivot point contained in faces 0 & 1
-        for(j = 0; j < (int)tempFaces[0].size(); j++)
+        if(this->CellZoneMesh != NULL)
           {
-          for(k = 0; k < (int)tempFaces[1].size(); k++)
+          for(size_t i = 0; i < this->CellZoneMesh->value.size(); i++)
             {
-            if(tempFaces[0][j] == tempFaces[1][k] && tempFaces[0][j] !=
-              firstFace[pointCount])
-              {
-              pivotPoint = tempFaces[0][j];
-              break;
-              }
+            this->CellZoneMesh->value[i]->SetPoints(points);
             }
           }
-        cellPoints.push_back(pivotPoint);
-        tempFaces[0].clear();
-        tempFaces[1].clear();
-        faceCount=0;
+        points->Delete();
         }
+      }
+    }
+  else
+    {
+    this->ClearMeshes();
 
-      //create the wedge cell and insert it into the mesh
-      vtkWedge * wedge = vtkWedge::New();
-      for(pCount = 0; pCount < (int)cellPoints.size(); pCount++)
-        {
-        wedge->GetPointIds()->SetId(pCount, cellPoints[pCount]);
-        }
-      cellZoneMesh->InsertNextCell(wedge->GetCellType(),
-                                   wedge->GetPointIds());
-      cellPoints.clear();
-      wedge->Delete();
-      firstFace.clear();
+    vtkStdString meshDir = *this->PathPrefix
+      + this->PolyMeshFacesDir->value[timeState] + "/polyMesh/";
+
+    //create paths to polyMesh files
+    vtkStdString facePath = meshDir + "faces";
+    vtkStdString ownerPath = meshDir + "owner";
+    vtkStdString neighborPath = meshDir + "neighbour";
+
+    //create the faces vector
+    intVectorVector* facePoints = this->ReadFacesFile(facePath.c_str());
+    if(facePoints == NULL)
+      {
+      return 0;
       }
+    this->UpdateProgress(0.2);
 
-    //OFpyramid | vtkPyramid
-    else if (totalPointCount == 16)
+    //read owner/neighbor and create the faces owner/facesOfCell vector
+    faceVectorVector* facesOfCell = new faceVectorVector;
+    if(!this->ReadOwnerNeighborFiles(facesOfCell, ownerPath.c_str(),
+      neighborPath.c_str()))
       {
-      foundDup = false;
+      delete facePoints;
+      delete facesOfCell;
+      return 0;
+      }
+    this->UpdateProgress(0.3);
 
-      //find quad
-      for(j = 0; j < (int)this->FacesOfCell->value[
-        cellZone[i]].size(); j++)  //each face
-        {
-        if((int)this->FacePoints->value[this->FacesOfCell->value[
-          cellZone[i]][j].faceIndex].size() == 4)
-          {
-          for(k = 0; k < (int)this->FacePoints->value[this->FacesOfCell->value[
-            cellZone[i]][j].faceIndex].size(); k++)
-            {
-            cellPoints.push_back(this->FacePoints->value[this->FacesOfCell->
-              value[cellZone[i]][j].faceIndex][k]);
-            }
-          break;
-          }
-        }
+    //get the points
+    vtkPoints* points = this->GetPoints(timeState);
+    if(points == NULL)
+      {
+      delete facePoints;
+      delete facesOfCell;
+      return 0;
+      }
+    this->UpdateProgress(0.4);
 
-      //compare first face points to second faces
-      for(j = 0; j < (int)cellPoints.size(); j++) //each point
+    // make internal mesh
+    // Create Internal Mesh only if required for display
+    if(this->GetPatchArrayStatus(GetPatchArrayName(0)))
+      {
+      this->InternalMesh = this->MakeInternalMesh(facesOfCell, facePoints,
+        points);
+      }
+    else
+      {
+      if(this->InternalMesh != NULL)
         {
-        for(k = 0; k < (int)this->FacePoints->value[this->FacesOfCell->value[
-          cellZone[i]][1].faceIndex].size(); k++)
-          {
-          if(cellPoints[j] == this->FacePoints->value[this->FacesOfCell->value[
-            cellZone[i]][1].faceIndex][k])
-            {
-            foundDup = true;
-            }
-          }
-        if(!foundDup)
-          {
-          cellPoints.push_back(this->FacePoints->value[this->FacesOfCell->value[
-            cellZone[i]][j].faceIndex][k]);
-          break;
-          }
+        this->InternalMesh->Delete();
+        this->InternalMesh = NULL;
         }
+      }
 
-      //create the pyramid cell and insert it into the mesh
-      vtkPyramid * pyramid = vtkPyramid::New();
-      for(pCount = 0; pCount < (int)cellPoints.size(); pCount++)
-        {
-        pyramid->GetPointIds()->SetId(pCount, cellPoints[pCount]);
-        }
-      cellZoneMesh->InsertNextCell(pyramid->GetCellType(),
-                                   pyramid->GetPointIds());
-      cellPoints.clear();
-      pyramid->Delete();
+    // read polyMesh/bondary
+    this->BoundaryDict = this->GatherBlocks("boundary", timeState, true);
+    if(this->BoundaryDict == NULL)
+      {
+      delete facePoints;
+      delete facesOfCell;
+      points->Delete();
+      vtkErrorMacro("Couldn't read polyMesh/boundary");
+      return 0;
       }
 
-    //OFtet | vtkTetrahedron
-    else if (totalPointCount == 12)
+    // create boundary mesh
+    this->BoundaryMesh = this->GetBoundaryMesh(this->BoundaryDict, facePoints,
+      points);
+    if(this->BoundaryMesh == NULL)
       {
-      foundDup = false;
+      delete facePoints;
+      delete facesOfCell;
+      points->Delete();
+      return 0;
+      }
 
-      //grab first face
-      for(j = 0; j < (int)this->FacePoints->value[this->FacesOfCell->value[
-        cellZone[i]][0].faceIndex].size(); j++)
+    // read and construct zones
+    if(ReadZones)
+      {
+      this->PointZoneMesh = new unstructuredGridVector;
+      if(!this->GetPointZoneMesh(this->PointZoneMesh, points, timeState))
+        {
+        for(size_t i = 0; i < this->PointZoneMesh->value.size(); i++)
+          {
+          this->PointZoneMesh->value[i]->Delete();
+          }
+        delete this->PointZoneMesh;
+        this->PointZoneMesh = NULL;
+        delete facePoints;
+        delete facesOfCell;
+        points->Delete();
+        return 0;
+        }
+      if(this->PointZoneMesh->value.size() == 0)
         {
-        cellPoints.push_back(this->FacePoints->value[this->FacesOfCell->value[
-          cellZone[i]][0].faceIndex][j]);
+        delete this->PointZoneMesh;
+        this->PointZoneMesh = NULL;
         }
 
-      //compare first face points to second faces
-      for(j = 0; j < (int)cellPoints.size(); j++) //each point
+      this->FaceZoneMesh = new unstructuredGridVector;
+      if(!this->GetFaceZoneMesh(this->FaceZoneMesh, facePoints, points,
+        timeState))
         {
-        for(k = 0; k < (int)this->FacePoints->value[this->FacesOfCell->value[
-          cellZone[i]][1].faceIndex].size(); k++)
+        for(size_t i = 0; i < this->FaceZoneMesh->value.size(); i++)
           {
-          if(cellPoints[j] == this->FacePoints->value[this->FacesOfCell->value[
-            cellZone[i]][1].faceIndex][k])
-            {
-            foundDup = true;
-            }
+          this->FaceZoneMesh->value[i]->Delete();
           }
-        if(!foundDup)
+        delete this->FaceZoneMesh;
+        this->FaceZoneMesh = NULL;
+        delete facePoints;
+        delete facesOfCell;
+        points->Delete();
+        return 0;
+        }
+      if(this->FaceZoneMesh->value.size() == 0)
+        {
+        delete this->FaceZoneMesh;
+        this->FaceZoneMesh = NULL;
+        }
+          
+      this->CellZoneMesh = new unstructuredGridVector;
+      if(!this->GetCellZoneMesh(this->CellZoneMesh, facesOfCell, facePoints,
+         points, timeState))
+        {
+        for(size_t i = 0; i < this->CellZoneMesh->value.size(); i++)
           {
-          cellPoints.push_back(this->FacePoints->value[this->FacesOfCell->value[
-            cellZone[i]][j].faceIndex][k]);
-          break;
+          this->CellZoneMesh->value[i]->Delete();
           }
+        delete this->CellZoneMesh;
+        this->CellZoneMesh = NULL;
+        delete facePoints;
+        delete facesOfCell;
+        points->Delete();
+        return 0;
         }
-
-      //create the wedge cell and insert it into the mesh
-      vtkTetra * tetra = vtkTetra::New();
-      for(pCount = 0; pCount < (int)cellPoints.size(); pCount++)
+      if(this->CellZoneMesh->value.size() == 0)
         {
-        tetra->GetPointIds()->SetId(pCount, cellPoints[pCount]);
+        delete this->CellZoneMesh;
+        this->CellZoneMesh = NULL;
         }
-      cellZoneMesh->InsertNextCell(tetra->GetCellType(),
-                                   tetra->GetPointIds());
-      cellPoints.clear();
-      tetra->Delete();
       }
 
-    //erronous cells
-    else if(totalPointCount == 0)
-      {
-      vtkWarningMacro("Warning: No points in cell.");
-      }
+    delete facesOfCell;
+    delete facePoints;
+    points->Delete();
+    }
+  this->UpdateProgress(0.5);
 
-    //OFpolyhedron || vtkConvexPointSet
-    else
+  if(updateVariables)
+    {
+    // read field data variables into Internal/Boundary meshes
+    for(int i = 0; i < (int)this->TimeStepData->value.size(); i++)
       {
-      vtkWarningMacro("Warning: Polyhedral Data is very Slow!");
-      foundDup = false;
-
-      //grab face 0
-      for(j = 0; j < (int)this->FacePoints->value[this->FacesOfCell->value[
-        cellZone[i]][0].faceIndex].size(); j++)
+      if(!this->GetVariableAtTimestep(this->InternalMesh, this->BoundaryMesh,
+        this->BoundaryDict, this->TimeStepData->value[i].c_str(), timeState))
         {
-        firstFace.push_back(this->FacePoints->value[
-          this->FacesOfCell->value[i][0].faceIndex][j]);
+        // error occurred while reading a field data variable, but continue
+        // to the next variable
+        continue;
         }
+      this->UpdateProgress(0.5
+        + 0.5 * ((float)(i + 1) / (float)this->TimeStepData->value.size()));
+      }
 
-      //ADD FIRST FACE TO CELL POINTS
-      for(j =0; j < (int)firstFace.size(); j++)
+    // read lagrangian mesh and fields
+    if(this->LagrangianTimeStepData->value.size() > 0
+      && this->GetPatchArrayStatus("Lagrangian Particles"))
+      {
+      // construct lagrangian mesh
+      lagrangianMesh = this->MakeLagrangianMesh(timeState);
+      if(lagrangianMesh == NULL)
         {
-        cellPoints.push_back(firstFace[j]);
+        return 0;
         }
-      //j = 1 skip firstFace
-      for(j = 1; j < (int) this->FacesOfCell->value[
-           cellZone[i]].size(); j++)
+
+      // read lagrangian variables
+      for(size_t i = 0; i < this->LagrangianTimeStepData->value.size(); i++)
         {
-        //remove duplicate points from faces
-        for(k = 0; k < (int)this->FacePoints->value[
-            this->FacesOfCell->value[i][j].faceIndex].size(); k++)
+        if(!this->GetLagrangianVariableAtTimestep(lagrangianMesh,
+          this->LagrangianTimeStepData->value[i].c_str(), timeState))
           {
-          for(l = 0; l < (int)cellPoints.size(); l++);
-            {
-            if(cellPoints[l] == this->FacePoints->value[this->FacesOfCell->
-                                value[cellZone[i]][j].faceIndex][k])
-              {
-              foundDup = true;
-              }
-            }
-          if(!foundDup)
-            {
-            cellPoints.push_back(this->FacePoints->value[this->FacesOfCell->
-                                 value[cellZone[i]][j].faceIndex][k]);
-            foundDup = false;
-            }
+          // error occurred while reading a field data variable, but continue
+          // to the next variable
+          continue;
           }
         }
+      }
+    }
 
-      //create the poly cell and insert it into the mesh
-      vtkConvexPointSet * poly = vtkConvexPointSet::New();
-      poly->GetPointIds()->SetNumberOfIds(cellPoints.size());
-      for(pCount = 0; pCount < (int)cellPoints.size(); pCount++)
-        {
-        poly->GetPointIds()->SetId(pCount, cellPoints[pCount]);
-        }
-      cellZoneMesh->InsertNextCell(poly->GetCellType(),
-                                   poly->GetPointIds());
-      cellPoints.clear();
-      firstFace.clear();
-      poly->Delete();
+  // assign group to each of data and point/face/cell-Zones
+  int nGroups = 1 + ((lagrangianMesh != NULL
+    && lagrangianMesh->GetPointData()->GetNumberOfArrays() > 0) ? 1 : 0)
+    + (this->PointZoneMesh != NULL ? 1 : 0)
+    + (this->FaceZoneMesh != NULL ? 1 : 0)
+    + (this->CellZoneMesh != NULL ? 1 : 0);
+  output->SetNumberOfGroups(nGroups);
+
+  // set internal mesh/data as output
+  //output->SetNumberOfDataSets(0,output->GetNumberOfDataSets(0));
+
+  // Add Internal Mesh to final output only if selected for display
+  if(GetPatchArrayStatus(GetPatchArrayName(0)))
+    {
+    output->SetDataSet(0, output->GetNumberOfDataSets(0), this->InternalMesh);
+    }
+
+  // set boundary meshes/data as output
+  for(size_t j = 0; j < this->BoundaryMesh->value.size(); j++)
+    {
+    output->SetDataSet(0, output->GetNumberOfDataSets(0),
+      this->BoundaryMesh->value[j]);
+    if(!this->CacheMesh)
+      {
+      this->BoundaryMesh->value[j]->Delete();
       }
     }
-  //set cell zone points
-  cellZoneMesh->SetPoints(Points);
-  input->close();
-  delete input;
-  vtkDebugMacro(<<"Cell zone mesh created");
-  return cellZoneMesh;
-}
 
-void vtkOpenFOAMReader::CreateDataSet(vtkMultiBlockDataSet *output)
-{
-  int timeState = this->TimeStep;
-  //create paths to polyMesh files
-  vtkstd::string boundaryPath = this->PathPrefix->value +
-                                this->PolyMeshFacesDir->value[timeState] +
-                                "/polyMesh/boundary";
-  vtkstd::string facePath = this->PathPrefix->value +
-                            this->PolyMeshFacesDir->value[timeState] +
-                            "/polyMesh/faces";
-  vtkstd::string ownerPath = this->PathPrefix->value +
-                             this->PolyMeshFacesDir->value[timeState] +
-                             "/polyMesh/owner";
-  vtkstd::string neighborPath = this->PathPrefix->value +
-                                this->PolyMeshFacesDir->value[timeState] +
-                                "/polyMesh/neighbour";
-  //create the faces vector
-  this->ReadFacesFile(facePath.c_str());
-
-  //create the faces owner vector
-  this->ReadOwnerFile(ownerPath.c_str());
-
-  //create the faces neighbor vector
-  this->ReadNeighborFile(neighborPath.c_str());
-
-  //create a vector containing a faces of each cell
-  this->CombineOwnerNeigbor();
-
-  this->GetPoints(timeState); //get the points
-
-  //get the names of the regions
-  //this->BoundaryNames->value =
-  //  this->GatherBlocks("boundary", timeState)->value;
-  //this->PointZoneNames->value =
-  //  this->GatherBlocks("pointZones", timeState)->value;
-  //this->FaceZoneNames->value =
-  //  this->GatherBlocks("faceZones", timeState)->value;
-  //this->CellZoneNames->value =
-  //  this->GatherBlocks("cellZones", timeState)->value;
-
-  //get the names of the regions
-  this->BoundaryNames =
-    this->GatherBlocks("boundary", timeState);
-  this->PointZoneNames =
-    this->GatherBlocks("pointZones", timeState);
-  this->FaceZoneNames =
-    this->GatherBlocks("faceZones", timeState);
-  this->CellZoneNames =
-    this->GatherBlocks("cellZones", timeState);
-
-  int numBoundaries = this->BoundaryNames->value.size();
-  int numPointZones = this->PointZoneNames->value.size();
-  int numFaceZones = this->FaceZoneNames->value.size();
-  int numCellZones = this->CellZoneNames->value.size();
-
-  //Internal Mesh
-  vtkUnstructuredGrid * internalMesh = this->MakeInternalMesh();
-  for(int i = 0; i < (int)this->TimeStepData->value.size(); i++)
-    {
-    vtkDoubleArray * data =
-      this->GetInternalVariableAtTimestep(this->TimeStepData->value[i].c_str(),
-                                          timeState);
-    if(data->GetSize() > 0)
+  if(!this->CacheMesh)
+    {
+    // First check if these meshes exist before deleting them....
+    if(this->BoundaryMesh != NULL)
       {
-      data->SetName(this->TimeStepData->value[i].c_str());
-      internalMesh->GetCellData()->AddArray(data);
+      delete this->BoundaryMesh;
+      this->BoundaryMesh = NULL;
       }
-    data->Delete();
+    if(this->InternalMesh != NULL)
+      {
+      this->InternalMesh->Delete();
+      this->InternalMesh = NULL;
+      }
+    delete this->BoundaryDict;
+    this->BoundaryDict = NULL;
+    this->FaceOwner->Delete();
+    this->FaceOwner = NULL;
     }
-  //output->SetNumberOfDataSets(0,output->GetNumberOfDataSets(0));
-  output->SetDataSet(0, output->GetNumberOfDataSets(0), internalMesh);
 
-  //Boundary Meshes
-  for(int i = 0; i < (int)numBoundaries; i++)
+  int groupI = 0;
+
+  // set lagrangian mesh as output
+  if(lagrangianMesh != NULL
+    && lagrangianMesh->GetPointData()->GetNumberOfArrays() > 0)
     {
-    vtkUnstructuredGrid * boundaryMesh = this->GetBoundaryMesh(timeState, i);
-    for(int j = 0; j < (int)this->TimeStepData->value.size(); j++)
+    groupI++;
+    output->SetDataSet(groupI, output->GetNumberOfDataSets(groupI),
+      lagrangianMesh);
+    lagrangianMesh->Delete();
+    }
+
+  // set Zone Meshes as output
+  if(this->PointZoneMesh != NULL)
+    {
+    groupI++;
+    for(size_t i = 0; i < this->PointZoneMesh->value.size(); i++)
       {
-      vtkDoubleArray * data =
-        this->GetBoundaryVariableAtTimestep(i, TimeStepData->value[j].c_str(),
-                                            timeState, internalMesh);
-      if(data->GetSize() > 0)
+      output->SetDataSet(groupI, output->GetNumberOfDataSets(groupI),
+        this->PointZoneMesh->value[i]);
+      if(!this->CacheMesh)
         {
-        data->SetName(this->TimeStepData->value[j].c_str());
-        boundaryMesh->GetCellData()->AddArray(data);
+        this->PointZoneMesh->value[i]->Delete();
         }
-      data->Delete();
       }
-    output->SetDataSet(0, output->GetNumberOfDataSets(0), boundaryMesh);
-    boundaryMesh->Delete();
+    if(!this->CacheMesh)
+      {
+      delete this->PointZoneMesh;
+      this->PointZoneMesh = NULL;
+      }
     }
 
-  internalMesh->Delete();
-  this->FaceOwner->Delete();
-  //Zone Meshes
-  for(int i = 0; i < (int)numPointZones; i++)
+  if(this->FaceZoneMesh != NULL)
     {
-    vtkUnstructuredGrid * pointMesh = this->GetPointZoneMesh(timeState, i);
-    output->SetDataSet(0, output->GetNumberOfDataSets(0),
-                       pointMesh);
-    pointMesh->Delete();
+    groupI++;
+    for(size_t i = 0; i < this->FaceZoneMesh->value.size(); i++)
+      {
+      output->SetDataSet(groupI, output->GetNumberOfDataSets(groupI),
+        this->FaceZoneMesh->value[i]);
+      if(!this->CacheMesh)
+        {
+        this->FaceZoneMesh->value[i]->Delete();
+        }
+      }
+    if(!this->CacheMesh)
+      {
+      delete this->FaceZoneMesh;
+      this->FaceZoneMesh = NULL;
+      }
     }
-  for(int i = 0; i < (int)numFaceZones; i++)
+
+  if(this->CellZoneMesh != NULL)
     {
-    vtkUnstructuredGrid * faceMesh = this->GetFaceZoneMesh(timeState, i);
-    output->SetDataSet(0, output->GetNumberOfDataSets(0),
-                       faceMesh);
-    faceMesh->Delete();
+    groupI++;
+    for(size_t i = 0; i < this->CellZoneMesh->value.size(); i++)
+      {
+      output->SetDataSet(groupI, output->GetNumberOfDataSets(groupI),
+        this->CellZoneMesh->value[i]);
+      if(!this->CacheMesh)
+        {
+        this->CellZoneMesh->value[i]->Delete();
+        }
+      }
+    if(!this->CacheMesh)
+      {
+      delete this->CellZoneMesh;
+      this->CellZoneMesh = NULL;
+      }
     }
-  for(int i = 0; i < (int)numCellZones; i++)
+
+  if(this->CacheMesh)
     {
-    vtkUnstructuredGrid * cellMesh = this->GetCellZoneMesh(timeState, i);
-    output->SetDataSet(0, output->GetNumberOfDataSets(0),
-                       cellMesh);
-    cellMesh->Delete();
+    this->OldTimeStep = timeState;
     }
-  return;
-}
-stdString * vtkOpenFOAMReader::GetLine(ifstream * file)
-{
-  char tempChar;
-  stdString * buffer = new stdString;
-  while(file->peek() != '\n')
+  else
     {
-    file->get(tempChar);
-    buffer->value += tempChar;
+    this->OldTimeStep = -1;
     }
-  file->get(tempChar); //throw out \n
-  return buffer;
+
+  return 1;
 }

--- ParaView3.orig/VTK/IO/vtkOpenFOAMReader.h	2007-11-07 21:51:56.000000000 +0100
+++ ParaView3/VTK/IO/vtkOpenFOAMReader.h	2007-10-22 13:22:01.000000000 +0200
@@ -37,16 +37,18 @@
 } face;
 
 class vtkUnstructuredGrid;
+class vtkPolyData;
 class vtkPoints;
 class vtkIntArray;
-class vtkFloatArray;
-class vtkDoubleArray;
+class vtkStringArray;
 class vtkDataArraySelection;
-struct stdString;
+class vtkCallbackCommand;
+class vtkStdString;
 struct stringVector;
 struct intVector;
 struct intVectorVector;
 struct faceVectorVector;
+struct unstructuredGridVector;
 
 class VTK_IO_EXPORT vtkOpenFOAMReader : public vtkMultiBlockDataSetAlgorithm
 {
@@ -93,16 +95,102 @@
   void DisableAllCellArrays();
   void EnableAllCellArrays();
 
+  // Description:
+  // Get the number of point arrays available in the input.
+  int GetNumberOfPointArrays(void);
+
+  // Description:
+  // Get/Set whether the point array with the given name is to
+  // be read.
+  int GetPointArrayStatus(const char* name);
+  void SetPointArrayStatus(const char* name, int status);
+
+  // Description:
+  // Get the name of the  point array with the given index in
+  // the input.
+  const char* GetPointArrayName(int index);
+
+  // Description:
+  // Turn on/off all point arrays.
+  void DisableAllPointArrays();
+  void EnableAllPointArrays();
+
+  // Description:
+  // Get the number of Patches (inlcuding Internal Mesh) available in the input.
+  int GetNumberOfPatchArrays(void);
+
+  // Description:
+  // Get/Set whether the Patch with the given name is to
+  // be read.
+  int GetPatchArrayStatus(const char* name);
+  void SetPatchArrayStatus(const char* name, int status);
+
+  // Description:
+  // Get the name of the Patch with the given index in
+  // the input.
+  const char* GetPatchArrayName(int index);
+
+  // Description:
+  // Turn on/off all Patches including the Internal Mesh.
+  void DisableAllPatchArrays();
+  void EnableAllPatchArrays();
+
+  // Selection observer callback functions
+  static void SelectionModifiedCallback
+  (
+    vtkObject* caller,
+    unsigned long eid,
+    void* clientdata,
+    void* calldata
+  );
+
+  void SelectionModified();
+
+  // Description:
+  // Set/Get whether mesh is to be cached.
+  vtkSetMacro(CacheMesh, int);
+  vtkGetMacro(CacheMesh, int);
+  vtkBooleanMacro(CacheMesh, int);
+
+  // Option to accumulate patches and fields across time steps
+  // Description:
+  // Set/Get whether the patches are listed even if size = 0
+  vtkSetMacro(KeepPatches,int);
+  vtkGetMacro(KeepPatches,int);
+  vtkBooleanMacro(KeepPatches, int);
+
+  // Description:
+  // Set/Get whether zones will be read.
+  vtkSetMacro(ReadZones, int);
+  vtkGetMacro(ReadZones, int);
+  vtkBooleanMacro(ReadZones, int);
 
 protected:
+  // for caching mesh
+  int CacheMesh;
+  // for reading point/face/cell-Zones
+  int ReadZones;
+
+  // for accumulating patches across time steps
+  int KeepPatches;
+
   vtkOpenFOAMReader();
   ~vtkOpenFOAMReader();
-  int RequestData(vtkInformation *,
-    vtkInformationVector **, vtkInformationVector *);
-  int RequestInformation(vtkInformation *,
-    vtkInformationVector **, vtkInformationVector *);
+  int RequestData(vtkInformation *, vtkInformationVector **,
+    vtkInformationVector *);
+  int RequestInformation(vtkInformation *, vtkInformationVector **,
+    vtkInformationVector *);
 
 private:
+
+  struct vtkFoamError;
+  struct vtkFoamToken;
+  struct vtkFoamFile;
+  struct vtkFoamIOobject;
+  struct vtkFoamEntryValue;
+  struct vtkFoamEntry;
+  struct vtkFoamDict;
+
   vtkOpenFOAMReader(const vtkOpenFOAMReader&);  // Not implemented.
   void operator=(const vtkOpenFOAMReader&);  // Not implemented.
 
@@ -113,52 +201,72 @@
   int TimeStep;
   int TimeStepRange[2];
   double * Steps;
-  bool RequestInformationFlag;
-  int StartFace;
 
-  stdString * Path;
-  stdString * PathPrefix;
-  stringVector * TimeStepData;
+  // Selection Observer for translating GUI Changes to Mesh Updates
+  vtkCallbackCommand * SelectionObserver;
+
+  // DataArraySelection for Patch / Region Data
+  vtkDataArraySelection * PatchDataArraySelection;
   vtkDataArraySelection * CellDataArraySelection;
-  intVectorVector * FacePoints;
-  intVectorVector * FacesOwnerCell;
-  intVectorVector * FacesNeighborCell;
-  faceVectorVector * FacesOfCell;
-  vtkPoints * Points;
+  vtkDataArraySelection * PointDataArraySelection;
+
+  // Previous and Current Status of each of the selection arrays (in
+  // Bitwise form). Maximum 31 patches or fields support as of now
+  // (since we use LSB for overflow handling)
+  unsigned long CellSelectionOldStatus;
+  unsigned long PointSelectionOldStatus;
+  unsigned long PatchSelectionOldStatus;
+  unsigned long CellSelectionStatus;
+  unsigned long PointSelectionStatus;
+  unsigned long PatchSelectionStatus;
+
+  vtkStdString* OldFileName;
+  vtkStdString* PathPrefix;
+  vtkStringArray *TimeNames;
+  stringVector * TimeStepData;
+  stringVector * LagrangianTimeStepData;
+
   vtkIdType NumCells;
-  vtkIdType NumFaces;
   vtkIntArray * FaceOwner;
-  //vtkIntArray * FaceNeighbor;
   stringVector * PolyMeshPointsDir;
   stringVector * PolyMeshFacesDir;
-  vtkIdType NumPoints;
-  intVector * SizeOfBoundary;
-  stringVector * BoundaryNames;
-  stringVector * PointZoneNames;
-  stringVector * FaceZoneNames;
-  stringVector * CellZoneNames;
-  int NumBlocks;
-
-  void CombineOwnerNeigbor();
-  vtkUnstructuredGrid * MakeInternalMesh();
-  double ControlDictDataParser(const char *);
-  void ReadControlDict ();
-  void GetPoints (int);
-  void ReadFacesFile (const char *);
-  void ReadOwnerFile(const char *);
-  void ReadNeighborFile(const char *);
+
+  // for caching mesh
+  int OldTimeStep;
+  vtkUnstructuredGrid* InternalMesh;
+  unstructuredGridVector* BoundaryMesh;
+  vtkFoamDict* BoundaryDict;
+  unstructuredGridVector* PointZoneMesh;
+  unstructuredGridVector* FaceZoneMesh;
+  unstructuredGridVector* CellZoneMesh;
+
+  // member functions
+  void ClearMeshes();
+  void MakeTimeStepData();
+  void InsertCellToGrid(vtkUnstructuredGrid *, int, faceVectorVector *,
+    intVectorVector*);
+  vtkUnstructuredGrid * MakeInternalMesh(faceVectorVector *, intVectorVector *,
+    vtkPoints *);
+  bool ListTimeDirectoriesByControlDict(vtkFoamDict* dict);
+  bool ListTimeDirectoriesByInstances();
+  bool ReadControlDict(const char *);
+  vtkPoints* GetPoints (int);
+  intVectorVector* ReadFacesFile (const char *);
+  bool ReadOwnerNeighborFiles(faceVectorVector *, const char *, const char *);
   void PopulatePolyMeshDirArrays();
-  const char * GetDataType(const char *, const char *);
-  vtkDoubleArray * GetInternalVariableAtTimestep(const char *, int);
-  vtkDoubleArray * GetBoundaryVariableAtTimestep(int, const char *, int,
-                                                 vtkUnstructuredGrid *);
-  stringVector *GatherBlocks(const char *, int);
-  vtkUnstructuredGrid * GetBoundaryMesh(int, int);
-  vtkUnstructuredGrid * GetPointZoneMesh(int, int);
-  vtkUnstructuredGrid * GetFaceZoneMesh(int, int);
-  vtkUnstructuredGrid * GetCellZoneMesh(int, int);
-  void CreateDataSet(vtkMultiBlockDataSet *);
-  stdString * GetLine(ifstream *);
+  bool GetVariableAtTimestep(vtkUnstructuredGrid*, unstructuredGridVector*,
+    vtkFoamDict*, const char*, int);
+  vtkFoamDict* GatherBlocks(const char *, int, bool);
+  unstructuredGridVector* GetBoundaryMesh(vtkFoamDict*, intVectorVector *,
+    vtkPoints *);
+  bool GetPointZoneMesh(unstructuredGridVector*, vtkPoints*, int);
+  bool GetFaceZoneMesh(unstructuredGridVector *, intVectorVector *,
+    vtkPoints *, int);
+  bool GetCellZoneMesh(unstructuredGridVector *, faceVectorVector *,
+    intVectorVector *, vtkPoints *, int);
+  vtkPolyData* MakeLagrangianMesh(int);
+  bool GetLagrangianVariableAtTimestep(vtkPolyData*, const char*, int);
+  int CreateDataSet(vtkMultiBlockDataSet *, int);
 };
 
 #endif