From 62b352a6fce24bad9e11613493d920b5c94de453 Mon Sep 17 00:00:00 2001
From: Stefen Guzik <guzik2@llnl.gov>
Date: Wed, 23 Jul 2008 02:05:07 +0000
Subject: [PATCH] CGNS I/O rewrite - alpha version

---
 Common/Context.h        |    4 +
 Common/CreateFile.cpp   |    4 +-
 Common/DefaultOptions.h |    2 +
 Common/Options.cpp      |   10 +
 Common/Options.h        |    1 +
 Geo/CGNSOptions.h       |   39 +
 Geo/CustomContainer.h   |  579 ++++++++++++
 Geo/GEntity.h           |    2 +-
 Geo/GModel.h            |    7 +-
 Geo/GModelIO_CGNS.cpp   | 1140 +++++++++++++++--------
 Geo/MEdge.h             |   15 +-
 Geo/MFace.h             |   21 +-
 Geo/MNeighbour.h        | 1941 ---------------------------------------
 Geo/MVertex.h           |    2 +-
 Geo/MZone.cpp           |  396 ++++++++
 Geo/MZone.h             |  446 +++++++++
 Geo/MZoneBoundary.cpp   |  902 ++++++++++++++++++
 Geo/MZoneBoundary.h     |  266 ++++++
 Geo/Makefile            |   66 +-
 Geo/SVector3.h          |   17 +-
 20 files changed, 3476 insertions(+), 2384 deletions(-)
 create mode 100644 Geo/CGNSOptions.h
 create mode 100644 Geo/CustomContainer.h
 delete mode 100644 Geo/MNeighbour.h
 create mode 100644 Geo/MZone.cpp
 create mode 100644 Geo/MZone.h
 create mode 100644 Geo/MZoneBoundary.cpp
 create mode 100644 Geo/MZoneBoundary.h

diff --git a/Common/Context.h b/Common/Context.h
index e5ac9094b8..7ac14b4d99 100644
--- a/Common/Context.h
+++ b/Common/Context.h
@@ -9,6 +9,8 @@
 #include <vector>
 #include <string>
 
+#include "CGNSOptions.h"
+
 // Interface-independent context 
 
 class Context_T {
@@ -169,6 +171,8 @@ class Context_T {
     int smooth_normals, reverse_all_normals;
     double angle_smooth_normals;
     double allow_swap_edge_angle;
+    int zone_definition;
+    CGNSOptions cgns_options;
   } mesh;
 
   // post processing options 
diff --git a/Common/CreateFile.cpp b/Common/CreateFile.cpp
index 3dc3dc60d7..231b7ec186 100644
--- a/Common/CreateFile.cpp
+++ b/Common/CreateFile.cpp
@@ -173,7 +173,9 @@ void CreateOutputFile(const char *filename, int format)
     break;
 
   case FORMAT_CGNS:
-    GModel::current()->writeCGNS(name, CTX.mesh.scaling_factor);
+    GModel::current()->writeCGNS(name, CTX.mesh.zone_definition,
+                                 CTX.mesh.cgns_options,
+                                 CTX.mesh.scaling_factor);
     break;
 
   case FORMAT_MED:
diff --git a/Common/DefaultOptions.h b/Common/DefaultOptions.h
index 2c51f1716c..1c82d20947 100644
--- a/Common/DefaultOptions.h
+++ b/Common/DefaultOptions.h
@@ -1058,6 +1058,8 @@ StringXNumber MeshOptions_Number[] = {
     "Display faces of volume mesh?" },
   { F|O, "VolumeNumbers" , opt_mesh_volumes_num , 0. , 
     "Display volume mesh element numbers?" },
+  { F|O, "ZoneDefinition" , opt_mesh_zone_definition , 2. , 
+    "Method for defining a zone (0=single zone, 1=by partition, 2=by physical)" },
 
   { 0, NULL , NULL , 0. , NULL }
 } ;
diff --git a/Common/Options.cpp b/Common/Options.cpp
index 06ad1a4da1..cb81369174 100644
--- a/Common/Options.cpp
+++ b/Common/Options.cpp
@@ -5232,6 +5232,16 @@ double opt_mesh_color_carousel(OPT_ARGS_NUM)
   return CTX.mesh.color_carousel;
 }
 
+double opt_mesh_zone_definition(OPT_ARGS_NUM)
+{
+  if(action & GMSH_SET){
+    CTX.mesh.zone_definition = (int)val;
+    if(CTX.mesh.zone_definition < 0 || CTX.mesh.zone_definition > 2)
+      CTX.mesh.zone_definition = 0;
+  }
+  return CTX.mesh.zone_definition;
+}
+
 double opt_mesh_nb_nodes(OPT_ARGS_NUM)
 {
   double s[50];
diff --git a/Common/Options.h b/Common/Options.h
index c79b37877c..82a114d7d4 100644
--- a/Common/Options.h
+++ b/Common/Options.h
@@ -489,6 +489,7 @@ double opt_mesh_draw_skin_only(OPT_ARGS_NUM);
 double opt_mesh_save_all(OPT_ARGS_NUM);
 double opt_mesh_save_groups_of_nodes(OPT_ARGS_NUM);
 double opt_mesh_color_carousel(OPT_ARGS_NUM);
+double opt_mesh_zone_definition(OPT_ARGS_NUM);
 double opt_mesh_nb_nodes(OPT_ARGS_NUM);
 double opt_mesh_nb_triangles(OPT_ARGS_NUM);
 double opt_mesh_nb_quadrangles(OPT_ARGS_NUM);
diff --git a/Geo/CGNSOptions.h b/Geo/CGNSOptions.h
new file mode 100644
index 0000000000..aeb79f4e7b
--- /dev/null
+++ b/Geo/CGNSOptions.h
@@ -0,0 +1,39 @@
+// Gmsh - Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
+//
+// See the LICENSE.txt file for license information. Please report all
+// bugs and problems to <gmsh@geuz.org>.
+//
+// CGNSOptions.h - Copyright (C) 2008 S. Guzik, C. Geuzaine, J.-F. Remacle
+#ifndef _CGNSOPTIONS_H_
+#define _CGNSOPTIONS_H_
+
+struct CGNSOptions
+{
+  // Types
+  enum CGNSConnectivityNodeType {
+    OneToOne,
+    Generalized
+  };
+
+  // Data
+  int vectorDim;                        // Number of dimensions in a vector
+                                        // (only relevant for a 2D mesh)
+  int gridConnectivityLocation;         // Location of connectivity
+                                        // Same as type GridLocation_t in
+                                        // 'cgnslib.h'
+  CGNSConnectivityNodeType connectivityNodeType;
+                                        // Type of CGNS connectivity description
+  std::string baseName;
+  std::string zoneName;
+
+  // Default values
+  CGNSOptions() :
+    vectorDim(2),
+    gridConnectivityLocation(2),  // Assumes Vertex = 2 in cgnslib.h
+    connectivityNodeType(Generalized),
+    baseName("Base_0"),
+    zoneName("Zone_&I0&")
+  { }
+};
+
+#endif
diff --git a/Geo/CustomContainer.h b/Geo/CustomContainer.h
new file mode 100644
index 0000000000..c1553069fe
--- /dev/null
+++ b/Geo/CustomContainer.h
@@ -0,0 +1,579 @@
+// Gmsh - Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
+//
+// See the LICENSE.txt file for license information. Please report all
+// bugs and problems to <gmsh@geuz.org>.
+//
+// CustomContainer.h - Copyright (C) 2008 S. Guzik, C. Geuzaine, J.-F. Remacle
+#ifndef _CUSTOMCONTAINER_H_
+#define _CUSTOMCONTAINER_H_
+
+#include <cstdlib>
+#include <cstring>
+
+
+/*******************************************************************************
+ *
+ * This lightweight pool class attempts to reduce memory usage by having the
+ * unused elements save links to each other using whatever memory they have
+ * available.
+ *
+ * The elements (type T) need to implement the member functions:
+ *
+ * void set_pool_prev(void *const)  - save an address
+ * void *get_pool_prev()            - return the saved address
+ *
+ * For example,
+ *
+ *----------*
+
+struct Box
+{
+   int i;
+   int j;
+   void set_pool_prev(void *const p)
+   {
+      *reinterpret_cast<void**>(&i) = p;
+   }
+   void *get_pool_prev()
+   {
+      return *reinterpret_cast<void**>(&i);
+   }
+};
+
+ *----------*
+ *
+ * Ensure that the memory used for each type T is at least 64 bits and note
+ * that this Pool does not do any construction or destruction of the elements.
+ *
+ ******************************************************************************/
+
+namespace CCon  // "Custom Container"
+{
+
+template <typename T>
+class Block;
+
+
+/*==============================================================================
+ *
+ * Class Pool
+ *
+ * Purpose
+ * =======
+ *
+ *   Interface to the pool
+ *
+ *============================================================================*/
+
+template <typename T>
+class Pool
+{
+
+ public:
+
+  // Constructor
+  Pool(const unsigned _blockSize = 128)
+    : tailBlock(0), tailElement(0), blockSize(_blockSize)
+  { }
+
+  // Destructor
+  ~Pool(){
+    while(tailBlock) {
+      Block<T> *const block = tailBlock;
+      tailBlock = block->prev;
+      delete block;
+    }
+  }
+
+  // Get an element
+  void *get()
+  {
+    if(!tailElement) create_block();
+    void *const rval = tailElement;
+    tailElement =  static_cast<T*>(tailElement->get_pool_prev());
+    return rval;
+  }
+
+  // Return an element
+  void remit(T *const elem)
+  {
+    elem->set_pool_prev(tailElement);
+    tailElement = elem;
+  }
+
+ private:
+
+  // Data
+  Block<T> *tailBlock;
+  T *tailElement;
+  unsigned blockSize;
+
+  // Create a new block
+  void create_block()
+  {
+    tailBlock = new Block<T>(tailBlock, blockSize);
+    const unsigned back = blockSize - 1;
+    tailBlock->array[back].set_pool_prev(tailElement);
+    tailElement = &tailBlock->array[back];
+    for(int n = back; n--; ) {
+      T *const prev = tailElement--;
+      tailElement->set_pool_prev(prev);
+    }
+  }
+
+  // Copy and assignment are not permitted
+  Pool(const Pool&);
+  Pool &operator=(const Pool&);
+
+};
+
+
+/*==============================================================================
+ *
+ * Class Block
+ *
+ * Purpose
+ * =======
+ *
+ *   Implements a block of uninitialized memory for type T
+ *
+ *============================================================================*/
+
+template <typename T>
+class Block
+{
+
+  // Data
+  Block *const prev;
+  T *array;
+
+  // Constructor
+  Block(Block *const _prev, const unsigned blockSize)
+    : prev(_prev)
+  {
+    array = static_cast<T*>(std::malloc(sizeof(T)*blockSize));
+  }
+
+  // Destructor
+  ~Block()
+  {
+    free(array);
+  }
+
+  friend class Pool<T>;
+
+};
+
+
+/*******************************************************************************
+ *
+ * This allocator-like class is used to implement many small vectors which are
+ * sized according to the number of boundary vertices commonly expected to
+ * connect to a vertex in 2D and 3D zones:
+ *   2 - expected faces for vertex on zone edge in 2D
+ *   6 - expected faces for vertex on zone face in 3D
+ *   8 - overflow faces for vertex on zone face in 3D
+ *  16 - approximate faces for vertices on zone edges and corners in 3D
+ * Memory for the above sizes are pooled for reuse.  Larger sizes are permitted
+ * and will be allocated using 'new' and 'delete'.
+ *
+ ******************************************************************************/
+
+/*==============================================================================
+ *
+ * Class FaceAllocator
+ *
+ * Purpose
+ * =======
+ *
+ *   Lightweight allocator for class FaceVector
+ *
+ * Notes
+ * =====
+ *
+ *   - This lightweight allocator expects only primitive types.  Constructors
+ *     and destructors are not called for type T.
+ *   - The memory is encapsulated in the Face* structures.  Handles to these
+ *     structures use the data array.  A pointer to the structure is regained
+ *     using the handle and an offset.  This is mostly cautionary as the offset
+ *     is normally expected to be zero.
+ *   - No per-object data is stored as per normal requirements for allocators.
+ *     Critical OpenMP sections are defined since multiple threads can access
+ *     the allocator.
+ *
+ *============================================================================*/
+
+template <typename T>
+class FaceAllocator
+{
+
+ public:
+
+  // Memory structures
+  struct Face2
+  {
+    T faces[2];
+    void set_pool_prev(void *const p)
+    {
+      *reinterpret_cast<void**>(faces) = p;
+    }
+    void *get_pool_prev()
+    {
+      return *reinterpret_cast<void**>(faces);
+    }
+    ptrdiff_t get_offset() const
+    {
+      return reinterpret_cast<const char*>(this) -
+        reinterpret_cast<const char*>(faces);
+    }
+  };
+  struct Face6
+  {
+    T faces[6];
+    void set_pool_prev(void *const p)
+    {
+      *reinterpret_cast<void**>(faces) = p;
+    }
+    void *get_pool_prev()
+    {
+      return *reinterpret_cast<void**>(faces);
+    }
+    ptrdiff_t get_offset() const
+    {
+      return reinterpret_cast<const char*>(this) -
+        reinterpret_cast<const char*>(faces);
+    }
+  };
+  struct Face8
+  {
+    T faces[8];
+    void set_pool_prev(void *const p) 
+    {
+      *reinterpret_cast<void**>(faces) = p;
+    }
+    void *get_pool_prev()
+    {
+      return *reinterpret_cast<void**>(faces);
+    }
+    ptrdiff_t get_offset() const
+    {
+      return reinterpret_cast<const char*>(this) -
+        reinterpret_cast<const char*>(faces);
+    }
+  };
+  struct Face16
+  {
+    T faces[16];
+    void set_pool_prev(void *const p) 
+    {
+      *reinterpret_cast<void**>(faces) = p;
+    }
+    void *get_pool_prev()
+    {
+      return *reinterpret_cast<void**>(faces);
+    }
+    ptrdiff_t get_offset() const
+    {
+      return reinterpret_cast<const char*>(this) -
+        reinterpret_cast<const char*>(faces);
+    }
+  };
+
+  // Set offsets
+  static void set_offsets()
+  {
+    Face2 f2;
+    offset2 = f2.get_offset();
+    Face6 f6;
+    offset6 = f6.get_offset();
+    Face2 f8;
+    offset8 = f8.get_offset();
+    Face16 f16;
+    offset16 = f16.get_offset();
+  }
+
+  // Allocate the array
+  void allocate(const unsigned short nCapacity, T *&faces)
+  {
+#pragma omp critical (FaceAllocator_allocate)
+    {
+      switch(nCapacity) {
+      case 0:
+        faces = 0;
+        break;
+      case 2:
+        {
+          Face2 *f2 = static_cast<Face2*>(face2Pool.get());
+          faces = f2->faces;
+        }
+        break;
+      case 6:
+        {
+          Face6 *f6 = static_cast<Face6*>(face6Pool.get());
+          faces = f6->faces;
+        }
+        break;
+      case 8:
+        {
+          Face8 *f8 = static_cast<Face8*>(face8Pool.get());
+          faces = f8->faces;
+        }
+        break;
+      case 16:
+        {
+          Face16 *f16 = static_cast<Face16*>(face16Pool.get());
+          faces = f16->faces;
+        }
+        break;
+      default:
+        {
+          faces = static_cast<T*>(std::malloc(sizeof(T)*nCapacity +
+                                              sizeof(void*)));
+        }
+        break;
+      }
+    }
+//omp: end critical
+  }
+
+  // Grow the array of faces by 2, 6, 8, and then 2*nUsed.  Pools are used for
+  // arrays with size 2, 6, and 8.
+  void grow(unsigned short &nCapacity, T *&faces)
+  {
+#pragma omp critical (FaceAllocator_grow)
+    {
+      switch(nCapacity) {
+      case 0:
+        {
+          Face2 *f2 = static_cast<Face2*>(face2Pool.get());
+          faces = f2->faces;
+          nCapacity = 2;
+        }
+        break;
+      case 2:
+        {
+          Face6 *f6 = static_cast<Face6*>(face6Pool.get());
+          std::memcpy(f6->faces, faces, 2*sizeof(T));
+          Face2 *f2 = reinterpret_cast<Face2*>(faces + offset2);
+          face2Pool.remit(f2);
+          faces = f6->faces;
+          nCapacity = 6;
+        }
+        break;
+      case 6:
+        {
+          Face8 *f8 = static_cast<Face8*>(face8Pool.get());
+          std::memcpy(f8->faces, faces, 6*sizeof(T));
+          Face6 *f6 = reinterpret_cast<Face6*>(faces + offset6);
+          face6Pool.remit(f6);
+          faces = f8->faces;
+          nCapacity = 8;
+        }
+        break;
+      case 8:
+        {
+          Face16 *f16 = static_cast<Face16*>(face16Pool.get());
+          std::memcpy(f16->faces, faces, 8*sizeof(T));
+          Face8 *f8 = reinterpret_cast<Face8*>(faces + offset8);
+          face8Pool.remit(f8);
+          faces = f16->faces;
+          nCapacity = 16;
+        }
+        break;
+      // Allocate outside pool for more than 16 faces
+      case 16:
+        {
+          Face16 *f16 = reinterpret_cast<Face16*>(faces + offset16);
+          faces = static_cast<T*>(std::malloc(sizeof(T)*32 + sizeof(void*)));
+          std::memcpy(faces, f16->faces, 16*sizeof(T));
+          face16Pool.remit(f16);
+          nCapacity = 32;
+        }
+        break;
+      default:
+        {
+          T *newFace = static_cast<T*>(std::malloc(sizeof(T)*2*nCapacity +
+                                                   sizeof(void*)));
+          std::memcpy(newFace, faces, nCapacity*sizeof(T));
+          std::free(faces);
+          faces = newFace;
+          nCapacity *= 2;
+        }
+        break;
+      }
+    }
+//omp: end critical
+  }
+
+  // Deallocate an array
+  void deallocate(unsigned short &nCapacity, T *const faces)
+  {
+#pragma omp critical (FaceAllocator_deallocate)
+    {
+      switch(nCapacity) {
+      case 0:
+        break;
+      case 2:
+        {
+          Face2 *const f2 = reinterpret_cast<Face2*>(faces + offset2);
+          face2Pool.remit(f2);
+        }
+        break;
+      case 6:
+        {
+          Face6 *const f6 = reinterpret_cast<Face6*>(faces + offset6);
+          face6Pool.remit(f6);
+        }
+        break;
+      case 8:
+        {
+          Face8 *const f8 = reinterpret_cast<Face8*>(faces + offset8);
+          face8Pool.remit(f8);
+        }
+        break;
+      case 16:
+        {
+          Face16 *const f16 = reinterpret_cast<Face16*>(faces + offset16);
+          face16Pool.remit(f16);
+        }
+        break;
+      default:
+        {
+          std::free(faces);
+        }
+        break;
+      }
+    }
+//omp: end critical
+    nCapacity = 0;
+  }
+
+ private:
+
+  // Data
+  static Pool<Face2> face2Pool;
+  static Pool<Face6> face6Pool;
+  static Pool<Face8> face8Pool;
+  static Pool<Face16> face16Pool;
+  static ptrdiff_t offset2;
+  static ptrdiff_t offset6;
+  static ptrdiff_t offset8;
+  static ptrdiff_t offset16;
+  
+};
+
+// Definitions for static data members of class FaceAllocator
+template <typename T>
+Pool<typename FaceAllocator<T>::Face2> FaceAllocator<T>::face2Pool;
+template <typename T>
+Pool<typename FaceAllocator<T>::Face6> FaceAllocator<T>::face6Pool;
+template <typename T>
+Pool<typename FaceAllocator<T>::Face8> FaceAllocator<T>::face8Pool;
+template <typename T>
+Pool<typename FaceAllocator<T>::Face16> FaceAllocator<T>::face16Pool;
+template <typename T>
+ptrdiff_t FaceAllocator<T>::offset2 = 0;
+template <typename T>
+ptrdiff_t FaceAllocator<T>::offset6 = 0;
+template <typename T>
+ptrdiff_t FaceAllocator<T>::offset8 = 0;
+template <typename T>
+ptrdiff_t FaceAllocator<T>::offset16 = 0;
+
+
+/*******************************************************************************
+ *
+ * This class is like a std::vector but optimized for small vectors containing
+ * elements with only primitive data.
+ *
+ ******************************************************************************/
+
+/*==============================================================================
+ *
+ * Class FaceVector
+ *
+ * Purpose
+ * =======
+ *
+ *   Lightweight small vectors
+ *
+ * Notes
+ * =====
+ *
+ *   - The only way to add elements is by 'push_back'
+ *   - Erasing may reorder the elements.
+ *   - T must only contain primitive types
+ *
+ *============================================================================*/
+
+template <typename T>
+class FaceVector : public FaceAllocator<T>
+{
+
+ public:
+
+  // Constructor
+  FaceVector() : _size(0), _capacity(0) { }
+  // Unlike std::vector, the following sets the capacity to 'n'.  The size is
+  // still 0.
+  FaceVector(const unsigned short n) : _size(0)
+  {
+      _capacity = valid_capacity(n);
+      allocate(_capacity, faces);
+  }
+
+  // Destructor
+  ~FaceVector() { deallocate(_capacity, faces); }
+
+  // Index the vector
+  const T &operator[](const int i) const { return faces[i]; }
+  T &operator[](const int i) { return faces[i]; }
+
+  // Add element to end
+  T &push_back(const T val)
+  {
+    if(_size == _capacity) grow(_capacity, faces);
+    return faces[_size++] = val;
+  }
+
+  // Just increment the size (push an empty element)
+  T &push_empty()
+  {
+    if(_size == _capacity) grow(_capacity, faces);
+    return faces[_size++];
+  }
+
+  // Erase an element
+  void erase(const int i)
+  {
+    faces[i] = faces[--_size];
+  }
+
+  // Vector size and capacity
+  unsigned size() const { return _size; }
+
+  unsigned capacity() const { return _capacity; }
+
+ private:
+
+  // Data
+  T *faces;
+  unsigned short _size;
+  unsigned short _capacity;
+
+  // Get a valid capacity size (returns n above 16)
+  unsigned valid_capacity(unsigned n) const
+  {
+    if(n == 0) return 0;
+    if(n <= 2) return 2;
+    if(n <= 6) return 6;
+    if(n <= 8) return 8;
+    if(n <= 16) return 16;
+    return n;
+  }
+
+};
+
+}  // End of namespace CCon
+
+#endif
diff --git a/Geo/GEntity.h b/Geo/GEntity.h
index 3596aa28a5..ed2e064d51 100644
--- a/Geo/GEntity.h
+++ b/Geo/GEntity.h
@@ -237,7 +237,7 @@ class GEntity {
   unsigned int getNumMeshVertices() { return mesh_vertices.size(); }
 
   // get the mesh vertex at the given index
-  MVertex *getMeshVertex(unsigned int index) { return mesh_vertices[index]; }
+//  MVertex *getMeshVertex(unsigned int index) { return mesh_vertices[index]; }
 };
 
 class GEntityLessThan {
diff --git a/Geo/GModel.h b/Geo/GModel.h
index 90cffe1ff6..0b324c9411 100644
--- a/Geo/GModel.h
+++ b/Geo/GModel.h
@@ -22,6 +22,7 @@ class GEO_Internals;
 class OCC_Internals;
 class smooth_normals;
 class FieldManager;
+class CGNSOptions;
 
 // A geometric model. The model is a "not yet" non-manifold B-Rep.
 class GModel
@@ -222,6 +223,9 @@ class GModel
   std::set<int> &getMeshPartitions() { return meshPartitions; }
   std::set<int> &recomputeMeshPartitions();
 
+  // Number of partitions
+  int getNumMeshPartitions() const { return meshPartitions.size(); }
+
   // Deletes all the partitions
   void deleteMeshPartitions();
 
@@ -292,7 +296,8 @@ class GModel
 
   // CFD General Notation System files
   int readCGNS(const std::string &name);
-  int writeCGNS(const std::string &name, double scalingFactor=1.0);
+  int writeCGNS(const std::string &name, const int zoneDefinition,
+                const CGNSOptions &options, double scalingFactor=1.0);
 
   // Med "Modele d'Echange de Donnees" file format (the static routine
   // is allowed to load multiple models/meshes)
diff --git a/Geo/GModelIO_CGNS.cpp b/Geo/GModelIO_CGNS.cpp
index 16ba3476d5..35a8e39510 100644
--- a/Geo/GModelIO_CGNS.cpp
+++ b/Geo/GModelIO_CGNS.cpp
@@ -2,87 +2,91 @@
 //
 // See the LICENSE.txt file for license information. Please report all
 // bugs and problems to <gmsh@geuz.org>.
-// Copyright (C) 2006 S. Guzik, C. Geuzaine, J.-F. Remacle
 //
-// This program is free software; you can redistribute it and/or modify
-// it under the terms of the GNU General Public License as published by
-// the Free Software Foundation; either version 2 of the License, or
-// (at your option) any later version.
-//
-// This program is distributed in the hope that it will be useful,
-// but WITHOUT ANY WARRANTY; without even the implied warranty of
-// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-// GNU General Public License for more details.
-//
-// You should have received a copy of the GNU General Public License
-// along with this program; if not, write to the Free Software
-// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
-// USA.
-// 
-// Please report all bugs and problems to <gmsh@geuz.org>.
-
-#include <string.h>
+// GModelIO_CGNS.cpp - Copyright (C) 2008 S. Guzik, C. Geuzaine, J.-F. Remacle
+
+#ifdef _OPENMP
+#include <omp.h>
+#else
+#define omp_get_num_threads() 1
+#define omp_get_thread_num() 0
+#define omp_lock_t int
+#define omp_init_lock(x)
+#define omp_set_lock(x)
+#define omp_unset_lock(x)
+#define omp_destroy_lock(x)
+#endif
+
+#include <cstring>
 #include <iostream>  // DBG
-#include <list>
 #include <map>
-#include <set>
 #include <string>
 #include <vector>
+#include <queue>
 
 #include "GModel.h"
 #include "Message.h"
-#include "MElement.h"
-#include "MNeighbour.h"
-#include "MVertex.h"
+#include "MZone.h"
+#include "MZoneBoundary.h"
+#include "CGNSOptions.h"
 
 #if defined(HAVE_LIBCGNS)
 
 #include <cgnslib.h>
 
+//--Error function for the CGNS library
+
 int cgnsErr(const int cgIndexFile = -1)
 {
-  Msg::Error(cg_get_error());
+  Message::Error(cg_get_error());
   if(cgIndexFile != -1)
-    if(cg_close(cgIndexFile)) Msg::Error("Even unable to close CGNS file");
+    if(cg_close(cgIndexFile)) Message::Error("Unable to close CGNS file");
   return 0;
 }
 
 
 /*==============================================================================
- * File scope structures used for CGNS I/O
+ * Required types
  *============================================================================*/
 
-//--Structure describing the zones sharing a vertex or element
-
-struct ZoneConn_t {
-  unsigned int indexZone;               // Index of the zone
-  unsigned int indexConn;               // Index of the vertex or element used
-                                        // to describe the connectivity
-  ZoneConn_t(const unsigned int iz, const unsigned int ic) :
-    indexZone(iz), indexConn(ic) { }
-};
+typedef std::map<int, std::vector<GEntity*> > PhysGroupMap;
+                                        // Type providing a vector of entities
+                                        // for a physical group index
 
 //--Class to make C style CGNS name strings act like C++ types
 
-class CgnsNameStr {
-  private:
+class CGNSNameStr
+{
+ private:
   char name[33];
-  public:
+ public:
   // Constructors
-  CgnsNameStr() { name[0] = '\0'; }
-  CgnsNameStr(const char *const cstr) {
+  CGNSNameStr() { name[0] = '\0'; }
+  CGNSNameStr(const char *const cstr)
+  {
     std::strncpy(name, cstr, 32);
     name[32] = '\0';
   }
-  CgnsNameStr(const CgnsNameStr& cgs) { std::strcpy(name, cgs.name); }
+  CGNSNameStr(std::string &s)
+  {
+    std::strncpy(name, s.c_str(), 32);
+    name[32] = '\0';
+  }
+  CGNSNameStr(const int d, const char *const fmt = "%d") 
+  {
+    std::sprintf(name, fmt, d);
+  }
+  CGNSNameStr(const CGNSNameStr& cgs) { std::strcpy(name, cgs.name); }
   // Assignment
-  CgnsNameStr& operator=(const CgnsNameStr& cgs) {
+  CGNSNameStr &operator=(const CGNSNameStr& cgs)
+  {
     if ( &cgs != this ) {
       std::strcpy(name, cgs.name);
     }
     return *this;
   }
-  CgnsNameStr& operator=(const char *const cstr) {
+  CGNSNameStr &operator=(const char *const cstr)
+  {
     std::strncpy(name, cstr, 32);
     name[32] = '\0';
     return *this;
@@ -92,8 +96,86 @@ class CgnsNameStr {
   const char *c_str() const { return name; }
 };
 
+//--Provides a CGNS element type for a MSH element type and indicates the order
+//--in which the CGNS elements are to be written:
+//    3D first-order elements
+//    3D second-order elements
+//    2D first-order elements
+//    2D second-order elements
+//    1D first-order elements
+//    1D second-order elements
+//    MSH_NUM_TYPE+1 is used to place non-cgns elements last.
+static const int msh2cgns[MSH_NUM_TYPE][2] = {
+  {BAR_2,          16},
+  {TRI_3,          11},
+  {QUAD_4,         12},
+  {TETRA_4,         1},
+  {HEXA_8,          4},
+  {PENTA_6,         3},
+  {PYRA_5,          2},
+  {BAR_3,          17},
+  {TRI_6,          13},
+  {QUAD_9,         15},
+  {TETRA_10,        5},
+  {HEXA_27,        10},
+  {PENTA_18,        8},
+  {PYRA_14,         6},
+  {-1, MSH_NUM_TYPE+1},  // MSH_PNT (NODE in CGNS but not used herein)
+  {QUAD_8,         14},
+  {HEXA_20,         9},
+  {PENTA_15,        7},
+  {-1, MSH_NUM_TYPE+1},  // MSH_PYR_13
+  {-1, MSH_NUM_TYPE+1},  // MSH_TRI_9
+  {-1, MSH_NUM_TYPE+1},  // MSH_TRI_10
+  {-1, MSH_NUM_TYPE+1},  // MSH_TRI_12
+  {-1, MSH_NUM_TYPE+1},  // MSH_TRI_15
+  {-1, MSH_NUM_TYPE+1},  // MSH_TRI_15I
+  {-1, MSH_NUM_TYPE+1},  // MSH_TRI_21
+  {-1, MSH_NUM_TYPE+1},  // MSH_LIN_4
+  {-1, MSH_NUM_TYPE+1},  // MSH_LIN_5
+  {-1, MSH_NUM_TYPE+1},  // MSH_LIN_6
+  {-1, MSH_NUM_TYPE+1},  // MSH_TET_20
+  {-1, MSH_NUM_TYPE+1},  // MSH_TET_35
+  {-1, MSH_NUM_TYPE+1},  // MSH_TET_56
+  {-1, MSH_NUM_TYPE+1},  // MSH_TET_34
+  {-1, MSH_NUM_TYPE+1}   // MSH_TET_52
+};
+
+//--This functor allows for sorting of the element types according to the
+//--"write-order" in array 'msh2cgns'
+
+struct ElemSortCGNSType
+{
+  bool operator()(const int t0, const int t1)
+  {
+    if(zoneElemConn[t0].numElem > 0 && zoneElemConn[t1].numElem > 0)
+      return msh2cgns[t0][1] < msh2cgns[t1][1];
+    else if(zoneElemConn[t0].numElem > 0)
+      return true;
+    else
+      return false;
+  }
+  ElemSortCGNSType(const ElementConnectivity *const _zoneElemConn)
+    : zoneElemConn(_zoneElemConn)
+  { }
+ private:
+  const ElementConnectivity *const zoneElemConn;
+};
+
 
 /*==============================================================================
+ * Forward declarations
+ *============================================================================*/
+
+template <unsigned DIM>
+int write_CGNS_zones(GModel &model, const int zoneDefinition,
+                     const CGNSOptions &options,
+                     const double scalingFactor, const int vectorDim,
+                     const PhysGroupMap &group,
+                     const int cgIndexFile, const int cgIndexBase);
+
+
+/*******************************************************************************
  *
  * Routine readCGNS
  *
@@ -107,7 +189,7 @@ class CgnsNameStr {
  *
  *   name               - (I) file name
  *
- *============================================================================*/
+ ******************************************************************************/
 
 int GModel::readCGNS(const std::string &name)
 {
@@ -115,58 +197,39 @@ int GModel::readCGNS(const std::string &name)
 }
 
 
-/*==============================================================================
+/*******************************************************************************
  *
  * Routine writeCGNS
  *
  * Purpose
  * =======
  *
- *   Writes out mesh in CGNS format
+ *   Writes out the mesh in CGNS format
  *
  * I/O
  * ===
  *
  *   name               - (I) file name
+ *   zoneDefinition     - (I) how to define a zone
+ *   options            - (I) options specific to CGNS
  *   scalingFactor      - (I) scaling for coordinates
  *
- *============================================================================*/
+ ******************************************************************************/
 
-int GModel::writeCGNS(const std::string &name, double scalingFactor)
+int GModel::writeCGNS(const std::string &name, const int zoneDefinition,
+                      const CGNSOptions &options, double scalingFactor)
 {
 
-//----- Type definitions -----//
-
-  typedef std::map<int, std::vector<GEntity*> >::iterator zone_iterator;
-  typedef std::vector<MTriangle*>::const_iterator const_triangle_iterator;
-  typedef std::vector<MQuadrangle*>::const_iterator const_quadrangle_iterator;
-  typedef std::list<MVertex*> BoundaryVertList;
-  typedef std::set<MVertex*> InteriorVertSet;
-  typedef std::list<MElement*> BoundaryElemList;
-  typedef std::set<MElement*> InteriorElemSet;
-
-//----- Parameters -----//
-
   enum {
-    vertex,
-    edge,
-    face,
-    region
+    vertex = 0,
+    edge   = 1,
+    face   = 2,
+    region = 3
   };
 
-  int meshDim = 0;                      // Dimension of the mesh
-  std::map<int, std::vector<GEntity*> > groups[4];
-                                        // vector of entities that belong to
+  PhysGroupMap groups[4];               // vector of entities that belong to
                                         // each physical group (in each
                                         // dimension)
-  std::multimap<MVertex*, ZoneConn_t> sharedVertex;
-                                        // map of shared vertices describing the
-                                        // zones that share it and the index of
-                                        // the vertex in each of the zones.
-//   std::multimap<MEdge, ZoneConn_t> sharedEdge;
-//   std::multimap<MFace, ZoneConn_t> sharedFace;
-  MNeighbour meshNeighbours;            // Database of neighbours (computed for
-                                        // each zone)
 
 //--Get a list of groups in each dimension and each of the entities in that
 //--group.  The groups will define the zones.  If no 2D or 3D groups exist, just
@@ -186,7 +249,7 @@ int GModel::writeCGNS(const std::string &name, double scalingFactor)
         if((*it)->triangles.size() + (*it)->quadrangles.size() > 0)
           groups[face][1].push_back(*it);
       if(!groups[face].size()) {
-        Msg::Error("No mesh elements were found");
+        Message::Error("No mesh elements were found");
         return 0;
       }
     }
@@ -195,8 +258,10 @@ int GModel::writeCGNS(const std::string &name, double scalingFactor)
 //--Open the file and get ready to write the zones
 
   // Will this be a 2D or 3D mesh?
+  int meshDim = 2;
+  int vectorDim = 3;
   if(groups[region].size()) meshDim = 3;
-  else meshDim = 2;
+  else vectorDim = options.vectorDim;
 
   // open the file
   int cgIndexFile;
@@ -204,377 +269,646 @@ int GModel::writeCGNS(const std::string &name, double scalingFactor)
 
   // write the base node
   int cgIndexBase;
-  if(cg_base_write(cgIndexFile, "Base_0000", 2, 2, &cgIndexBase))
+  if(cg_base_write(cgIndexFile, "Base_0000", meshDim, meshDim, &cgIndexBase))
     return cgnsErr();
 
   // write information about who generated the mesh
   if(cg_goto(cgIndexFile, cgIndexBase, "end")) return(cgnsErr());
   if(cg_descriptor_write("About", "Created by Gmsh")) return cgnsErr();
 
-  unsigned int indexZone = 0;           // Index of a zone
-  const unsigned int nZone = groups[face].size() + groups[region].size();
-                                        // Total number of zones
-  std::vector<CgnsNameStr> zoneName(nZone);
-                                        // zone names
-  std::vector<BoundaryVertList> zoneBoundaryVert(nZone);
+  switch(meshDim) {
+  case 2:
+    write_CGNS_zones<2>(*this, zoneDefinition, options, scalingFactor,
+                        vectorDim, groups[face], cgIndexFile, cgIndexBase);
+    break;
+  case 3:
+    write_CGNS_zones<3>(*this, zoneDefinition, options, scalingFactor,
+                        vectorDim, groups[region], cgIndexFile, cgIndexBase);
+    break;
+  }
 
-/*--------------------------------------------------------------------*
- * Write the 3D zone information
- *--------------------------------------------------------------------*/
+//--Close the file
 
-// TODO: add
+  if(cg_close(cgIndexFile)) return cgnsErr();
 
-/*--------------------------------------------------------------------*
- * Write the 2D zone information
- *--------------------------------------------------------------------*/
+  return 0;
+}
 
-  for(zone_iterator itZone = groups[face].begin(); itZone != groups[face].end();
-      ++itZone) {
-    const int nFace = itZone->second.size();
-    std::cout << "Working on zone " << indexZone+1 << "  " << itZone->first
-              << std::endl;
-    std::sprintf(zoneName[indexZone].c_str(), "Zone_%d\0", indexZone+1);
-
-//--Compute the mesh neighbours for this group
-
-    meshNeighbours.add_elements_in_entities<GFace>(itZone->second.begin(),
-                                                   itZone->second.end());
-
-//--Load interior lists with all vertices and all elements
-
-    unsigned int numTri3 = 0;
-    unsigned int numTri6 = 0;
-    unsigned int numQuad4 = 0;
-    unsigned int numQuad8 = 0;
-    unsigned int numQuad9 = 0;
-
-    InteriorVertSet interiorVert;
-    BoundaryElemList boundaryElem;
-    InteriorElemSet interiorElem;
-
-    for(unsigned int iFace = 0; iFace != nFace; ++iFace) {
-      GFace* face = static_cast<GFace*>(itZone->second[iFace]);
-      const std::vector<MTriangle*>::const_iterator itEndTri =
-        face->triangles.end();
-      for(std::vector<MTriangle*>::const_iterator itTri =
-            face->triangles.begin(); itTri != itEndTri ; ++itTri) {
-        interiorElem.insert(*itTri);
-        interiorVert.insert((*itTri)->getVertex(0));
-        interiorVert.insert((*itTri)->getVertex(1));
-        interiorVert.insert((*itTri)->getVertex(2));
-        switch ( (*itTri)->getTypeForMSH() ) {
-        case MSH_TRI_3:
-          ++numTri3;
-          break;
-        case MSH_TRI_6:
-          ++numTri6;
-          interiorVert.insert((*itTri)->getVertex(3));
-          interiorVert.insert((*itTri)->getVertex(4));
-          interiorVert.insert((*itTri)->getVertex(5));
-          break;
-        }
+
+/*******************************************************************************
+ *
+ * Routine tranlateElementMSH2CGNS
+ *
+ * Purpose
+ * =======
+ *
+ *   Rewrites the mesh connectivity using CGNS node ordering
+ *
+ * I/O
+ * ===
+ *
+ *   zoneElemConn       - (I) connectivity for the zone using gmsh node ordering
+ *                            (see MElement.h)
+ *                        (O) connectivity using CGNS node ordering
+ *   
+ ******************************************************************************/
+
+void translateElementMSH2CGNS(ElementConnectivity *const zoneElemConn)
+{
+  const int maxVPE = 27;
+  // Location to write a MSH_TET_10 vertex to get a CGNS TETRA_10
+  static const int trTET10[10] = {
+    0, 1, 2, 3, 4, 5, 6, 7, 9, 8
+  };
+  // Location to write a MSH_HEX_20 vertex to get a CGNS HEXA_20
+  static const int trHEX20[20] = {
+    0, 1, 2, 3, 4, 5, 6, 7, 8, 11, 12, 9, 13, 10, 14, 15, 16, 19, 17, 18
+  };
+  // Location to write a MSH_HEX_27 vertex to get a CGNS HEXA_27
+  static const int trHEX27[27] = {
+    0, 1, 2, 3, 4, 5, 6, 7, 8, 11, 12, 9, 13, 10, 14, 15, 16, 19, 17, 18,
+    20, 21, 24, 22, 23, 25, 26
+  };
+  // Location to write a MSH_PRI_15 vertex to get a CGNS PENTA_15
+  static const int trPRI15[15] = {
+    0, 1, 2, 3, 4, 5, 6, 8, 9, 7, 10, 11, 12, 14, 13
+  };
+  // Location to write a MSH_PRI_18 vertex to get a CGNS PENTA_18
+  static const int trPRI18[18] = {
+    0, 1, 2, 3, 4, 5, 6, 8, 9, 7, 10, 11, 12, 14, 13,
+    15, 17, 16
+  };
+  // Location to write a MSH_PYR_13 vertex to get a CGNS PYRA_13
+  // PYRA_13 are not in the CGNS standard but if ever added they would have the
+  // same ordering as PYRA_14
+  static const int trPYR13[13] = {
+    0, 1, 2, 3, 4, 5, 8, 9, 6, 10, 7, 11, 12
+  };
+  // Location to write a MSH_PYR_14 vertex to get a CGNS PYRA_14
+  static const int trPYR14[14] = {
+    0, 1, 2, 3, 4, 5, 8, 9, 6, 10, 7, 11, 12,
+    13
+  };
+
+  int tmp[maxVPE];
+  for(int iElemType = 0; iElemType != MSH_NUM_TYPE; ++iElemType) {
+    int numVPE;
+    const int *tr;
+    if(zoneElemConn[iElemType].numElem > 0) {
+      switch(iElemType + 1) {
+      case MSH_TET_10:
+          numVPE = 10;
+          tr = trTET10;
+        break;
+      case MSH_HEX_20:
+          numVPE = 20;
+          tr = trHEX20;
+        break;
+      case MSH_HEX_27:
+          numVPE = 27;
+          tr = trHEX27;
+        break;
+      case MSH_PRI_15:
+          numVPE = 15;
+          tr = trPRI15;
+        break;
+      case MSH_PRI_18:
+          numVPE = 18;
+          tr = trPRI18;
+        break;
+      case MSH_PYR_13:  // Although not in the CGNS standard, I assume it would
+                        // be the same as PYRA_14 if ever added
+          numVPE = 13;
+          tr = trPYR13;
+        break;
+      case MSH_PYR_14:
+          numVPE = 14;
+          tr = trPYR14;
+        break;
+      default:
+        numVPE = 0;
+        tr = 0;
+        break;
       }
-      const std::vector<MQuadrangle*>::const_iterator itEndQuad =
-        face->quadrangles.end();
-      for(std::vector<MQuadrangle*>::const_iterator itQuad =
-            face->quadrangles.begin(); itQuad != itEndQuad ; ++itQuad) {
-        interiorElem.insert(*itQuad);
-        interiorVert.insert((*itQuad)->getVertex(0));
-        interiorVert.insert((*itQuad)->getVertex(1));
-        interiorVert.insert((*itQuad)->getVertex(2));
-        interiorVert.insert((*itQuad)->getVertex(3));
-        switch ( (*itQuad)->getTypeForMSH() ) {
-        case MSH_QUA_4:
-          ++numQuad4;
-          break;
-        case MSH_QUA_8:
-          ++numQuad8;
-          interiorVert.insert((*itQuad)->getVertex(4));
-          interiorVert.insert((*itQuad)->getVertex(5));
-          interiorVert.insert((*itQuad)->getVertex(6));
-          interiorVert.insert((*itQuad)->getVertex(7));
-          break;
-        case MSH_QUA_9:
-          ++numQuad9;
-          interiorVert.insert((*itQuad)->getVertex(4));
-          interiorVert.insert((*itQuad)->getVertex(5));
-          interiorVert.insert((*itQuad)->getVertex(6));
-          interiorVert.insert((*itQuad)->getVertex(7));
-          interiorVert.insert((*itQuad)->getVertex(8));
-          break;
+      if(numVPE > 0) {
+        int *aConn = &zoneElemConn[iElemType].connectivity[0];
+        for(int n = zoneElemConn[iElemType].numElem; n--; ) {
+          std::memcpy(tmp, aConn, numVPE*sizeof(int));
+          for(int iV = 0; iV != numVPE; ++iV) aConn[tr[iV]] = tmp[iV];
+          aConn += numVPE;
         }
       }
     }
+  }
+}
 
-//--Check edges for zone boundaries and move vertices and elements from interior
-//--sets to boundary lists
-
-    const MNeighbour::Edge_const_iterator itEndEdge =
-      meshNeighbours.edge_end();
-    for(MNeighbour::Edge_const_iterator itEdge = meshNeighbours.edge_begin();
-        itEdge != itEndEdge; ++itEdge) {
-      if(itEdge.num_neighbours() == 1) {  // Boundary edge
-        // Find the boundary element
-        MElement *element;
-        meshNeighbours.edge_neighbours(itEdge, &element);
-        // Move primary vertices from interior to boundary
-        InteriorVertSet::iterator itVert =
-          interiorVert.find(itEdge->getVertex(0));
-        if(itVert != interiorVert.end()) {
-          zoneBoundaryVert[indexZone].push_back(*itVert);
-          interiorVert.erase(itVert);
-        }
-        itVert = interiorVert.find(itEdge->getVertex(1));
-        if(itVert != interiorVert.end()) {
-          zoneBoundaryVert[indexZone].push_back(*itVert);
-          interiorVert.erase(itVert);
-        }
-        // If this is a higher-order element, find the edge on the boundary
-        // and then the mid-edge vertex
-        if(element->getNumVertices() != element->getNumPrimaryVertices()) {
-          // Second-order element
-          int iEPE = 0;
-          while(element->getEdge(iEPE) != *itEdge) ++iEPE;
-          // Add mid-vertex on edge iEPE
-          itVert = interiorVert.find
-            (element->getVertex(iEPE + element->getNumPrimaryVertices()));
-          if(itVert != interiorVert.end()) {
-            zoneBoundaryVert[indexZone].push_back(*itVert);
-            interiorVert.erase(itVert);
-          }
-        }
-      }
-    }
 
-//--Loop through all the boundary vertices
-
-    int numVert = 1;
-    BoundaryVertList::const_iterator itEndBVert =
-      zoneBoundaryVert[indexZone].end();
-    {
-      MElement **boElem = new MElement*[meshNeighbours.max_vertex_neighbours()];
-      for(BoundaryVertList::iterator itVert =
-            zoneBoundaryVert[indexZone].begin(); itVert != itEndBVert;
-          ++itVert) {
-        // Any elements with a vertex in the boundary vertex list should be
-        // moved to the boundary element list.
-        const int nElem = meshNeighbours.vertex_neighbours(*itVert, boElem);
-        for(int iElem = 0; iElem != nElem; ++iElem) {
-          InteriorElemSet::iterator itElem = interiorElem.find(boElem[iElem]);
-          if(itElem != interiorElem.end()) {
-            boundaryElem.push_back(boElem[iElem]);
-            interiorElem.erase(itElem);
-          }
-        }
-        // Renumber the boundary vertices
-        (*itVert)->setNum(numVert);
-        // Add all boundary vertices to the map of shared vertices
-        sharedVertex.insert(std::make_pair(*itVert, ZoneConn_t(indexZone,
-                                                               numVert)));
-        ++numVert;
-      }
-      delete[] boElem;
-    }
+/*******************************************************************************
+ *
+ * Routine get_zone_definition
+ *
+ * Purpose
+ * =======
+ *
+ *   Defines the next zone based on physicals or partitions and provides a name
+ *   for the zone
+ *
+ * I/O
+ * ===
+ *
+ *   model              - (I) gmsh model
+ *   zoneDefinition     - (I) how to define the zone (see enum in code)
+ *   options            - (I) options for CGNS I/O
+ *   numPartitions      - (I)
+ *   group              - (I) the group of physicals used to define the mesh
+ *   globalZoneIndex    - (I/O)
+ *   globalParittion    - (I/O)
+ *   globalItPhysical   - (I/O) a global scope iterator to the physicals
+ *   zoneIndex          - (O) index of the zone
+ *   partition          - (O) partition of the zone
+ *   itPhysicalBegin    - (O) begin physical for defining the zone
+ *   itPhysicalEnd      - (O) end physical for defining the zone
+ *   zoneName           - (O) name of the zone
+ *
+ ******************************************************************************/
+
+int get_zone_definition(GModel &model, const int zoneDefinition,
+                        const CGNSOptions &options,
+                        const int numPartitions, const PhysGroupMap &group,
+                        int &globalZoneIndex, int &globalPartition,
+                        PhysGroupMap::const_iterator &globalItPhysical,
+                        int &zoneIndex, int &partition,
+                        PhysGroupMap::const_iterator &itPhysicalBegin,
+                        PhysGroupMap::const_iterator &itPhysicalEnd,
+                        CGNSNameStr &zoneName)
+{
+  enum {
+    DefineZoneSingle = 0,
+    DefineZoneByPartition = 1,
+    DefineZoneByPhysical = 2
+  };
 
-//--Renumber the interior mesh vertices
-
-    InteriorVertSet::const_iterator itEndIVert = interiorVert.end();
-    for(InteriorVertSet::iterator itVert = interiorVert.begin();
-        itVert != itEndIVert; ++itVert) (*itVert)->setNum(numVert++);
-
-//--Write the zone node
-
-    int cgIndexZone;
-    int cgZoneSize[3];
-    cgZoneSize[0] = zoneBoundaryVert[indexZone].size() + interiorVert.size();
-    cgZoneSize[1] = boundaryElem.size() + interiorElem.size();
-    cgZoneSize[2] = zoneBoundaryVert[indexZone].size();
-    if(cg_zone_write(cgIndexFile, cgIndexBase, zoneName[indexZone].c_str(),
-                     cgZoneSize, Unstructured, &cgIndexZone)) return cgnsErr();
-
-//--Write the grid node
-
-    int cgIndexGrid;
-    if(cg_grid_write(cgIndexFile, cgIndexBase, cgIndexZone, "GridCoordinates",
-                     &cgIndexGrid)) return cgnsErr();
-
-//--Write the grid coordinates
-    
-    {
-      std::vector<double> coordx(cgZoneSize[0]);
-      std::vector<double> coordy(cgZoneSize[1]);
-      int iVert = 0;
-      for(BoundaryVertList::const_iterator itVert =
-            zoneBoundaryVert[indexZone].begin(); itVert != itEndBVert;
-          ++itVert) {
-        coordx[iVert] = (*itVert)->x()*scalingFactor;
-        coordy[iVert] = (*itVert)->y()*scalingFactor;
-        ++iVert;
-      }
-      for(InteriorVertSet::const_iterator itVert = interiorVert.begin();
-          itVert != itEndIVert; ++itVert) {
-        coordx[iVert] = (*itVert)->x()*scalingFactor;
-        coordy[iVert] = (*itVert)->y()*scalingFactor;
-        ++iVert;
-      }      
-      int cgIndexCoord;
-      if(cg_coord_write(cgIndexFile, cgIndexBase, cgIndexZone, RealDouble,
-                        "CoordinateX", &coordx[0], &cgIndexCoord))
-        return cgnsErr();
-      if(cg_coord_write(cgIndexFile, cgIndexBase, cgIndexZone, RealDouble,
-                        "CoordinateY", &coordy[0], &cgIndexCoord))
-        return cgnsErr();
-    }
+  int status = 0;
+  const char *_zoneName = "Partition";
 
-//--Write out the element connectivity
-    
-    {
-      int mixedElem = -1;
-      int numVPE;
-      int numElemData = 0;
-      ElementType_t cgElemType;
-      if(numTri3) {
-        ++mixedElem;
-        numVPE = 3;
-        numElemData += numTri3*3;
-        cgElemType = TRI_3;
-      }
-      if(numTri6) {
-        ++mixedElem;
-        numVPE = 6;
-        numElemData += numTri6*6;
-        cgElemType = TRI_6;
-      }
-      if(numQuad4) {
-        ++mixedElem;
-        numVPE = 4;
-        numElemData += numQuad4*4;
-        cgElemType = QUAD_4;
+//--Get indices for the zone
+
+#pragma omp critical (get_zone_definition)
+  {
+    switch(zoneDefinition) {
+    case DefineZoneSingle:
+      if(globalZoneIndex > 0) status = 1;  // No more zones
+      else {
+        partition = -1;
+        itPhysicalBegin = group.begin();
+        itPhysicalEnd = group.end();
       }
-      if(numQuad8) {
-        ++mixedElem;
-        numVPE = 8;
-        numElemData += numQuad8*8;
-        cgElemType = QUAD_8;
+      break;
+    case DefineZoneByPartition:
+      if(globalPartition > numPartitions) status = 1;   // No more zones
+      else {
+        partition = globalPartition++;
+        itPhysicalBegin = group.begin();
+        itPhysicalEnd = group.end();
       }
-      if(numQuad9) {
-        ++mixedElem;
-        numVPE = 9;
-        numElemData += numQuad9*9;
-        cgElemType = QUAD_9;
+      break;
+    case DefineZoneByPhysical:
+      if(globalItPhysical == group.end()) status = 1;  // No more zones
+      else {
+        partition = -1;
+        _zoneName = model.getPhysicalName(globalItPhysical->first).c_str();
+        itPhysicalBegin = globalItPhysical++;
+        itPhysicalEnd = globalItPhysical;
       }
-      if(mixedElem) {
-        numVPE = 0;
-        numElemData += cgZoneSize[1];
-        cgElemType = MIXED;
+      break;
+    }
+    zoneIndex = globalZoneIndex++;
+  }
+//omp: end critical
+
+//--Set the name for the zone
+
+  if(status == 0) {
+    std::string s = options.zoneName;
+    std::string::size_type r1 = s.find('&');
+    // Look for and replace "&&" sub-strings
+    while(r1 != std::string::npos) {
+      const std::string::size_type r2 = s.find('&', r1+1);
+      if(r2 == std::string::npos) {
+        s.replace(r1, s.length()-r1, "");
       }
-
-      std::vector<int> elementConn(numElemData);
-      int iElemData = 0;
-
-      // Load boundary elements
-      BoundaryElemList::const_iterator itEndBElem = boundaryElem.end();
-      if(mixedElem) {
-        for(BoundaryElemList::const_iterator itElem = boundaryElem.begin();
-            itElem != itEndBElem; ++itElem) {
-          switch((*itElem)->getTypeForMSH()) {
-          case MSH_TRI_3:
-            elementConn[iElemData++] = TRI_3;
-            break;
-          case MSH_TRI_6:
-            elementConn[iElemData++] = TRI_6;
-            break;
-          case MSH_QUA_4:
-            elementConn[iElemData++] = QUAD_4;
-            break;
-          case MSH_QUA_8:
-            elementConn[iElemData++] = QUAD_8;
-            break;
-          case MSH_QUA_9:
-            elementConn[iElemData++] = QUAD_9;
+      else {
+        const int rn = r2 - r1 + 1;
+        switch(s[r1+1]) {
+        case 'I':
+          // &I[0][%width]&
+          {
+            int iplus = 1;
+            if(s[r1+2] == '0') iplus = 0;
+            char fmt[6] = "%d";
+            const std::string::size_type f = s.find('%', r1+1);
+            if(f != std::string::npos && f < r2) {
+              int w = std::atoi(s.substr(f+1, r2-f-1).c_str());
+              if(w > 0 && w < 33) std::sprintf(fmt, "%%0%dd", w);
+            }
+            s.replace(r1, rn, CGNSNameStr(zoneIndex+iplus, fmt).c_str());
             break;
           }
-          const int nVPE = (*itElem)->getNumVertices();
-          for(int iVPE = 0; iVPE != nVPE; ++iVPE)
-            elementConn[iElemData++] = (*itElem)->getVertex(iVPE)->getNum();
+        case 'N':
+          // &N&
+          s.replace(r1, rn, _zoneName);
+          break;
+        default:
+          s.replace(r1, rn, "");
         }
       }
-      else {
-        const int nVPE = numVPE;
-        for(BoundaryElemList::const_iterator itElem = boundaryElem.begin();
-            itElem != itEndBElem; ++itElem)
-          for(int iVPE = 0; iVPE != nVPE; ++iVPE)
-            elementConn[iElemData++] = (*itElem)->getVertex(iVPE)->getNum();
-      }
+      if(s.length() > 1024) s = "";  // idiot/recursive check
+      r1 = s.find('&');
+    }
+    if(s.length() == 0) {  // If empty
+      s = "Zone_";
+      s += CGNSNameStr(zoneIndex+1).c_str();
+    }
+    zoneName = s.c_str();
+  }
 
-      // Load interior elements
-      InteriorElemSet::const_iterator itEndIElem = interiorElem.end();
-      if(mixedElem) {
-        for(InteriorElemSet::const_iterator itElem = interiorElem.begin();
-            itElem != itEndIElem; ++itElem) {
-          switch((*itElem)->getTypeForMSH()) {
-          case MSH_TRI_3:
-            elementConn[iElemData++] = TRI_3;
-            break;
-          case MSH_TRI_6:
-            elementConn[iElemData++] = TRI_6;
-            break;
-          case MSH_QUA_4:
-            elementConn[iElemData++] = QUAD_4;
-            break;
-          case MSH_QUA_8:
-            elementConn[iElemData++] = QUAD_8;
-            break;
-          case MSH_QUA_9:
-            elementConn[iElemData++] = QUAD_9;
-            break;
+  return status;
+}
+
+
+/*******************************************************************************
+ *
+ * Routine write_CGNS_zone
+ *
+ * Purpose
+ * =======
+ *
+ *   Writes the CGNS zones, splitting up the work between several processors if
+ *   threads are enabled.
+ *
+ * I/O
+ * ===
+ *
+ *   model              - (I) 
+ *   zoneDefinition     - (I) how to define a zone
+ *   options            - (I) CGNS specific options
+ *   scalingFactor
+ *   vectorDim          - (I) Dimensions of a vector (must be 3 for a 3D mesh,
+ *                            may be 2 or 3 for a 2D mesh)
+ *   group              - (I) Group of physicals and associated entities
+ *   cgIndexFile        - (I) index of the CGNS file node
+ *   cgIndexBase        - (I) index of the CGNS base node
+ *
+ ******************************************************************************/
+
+/*--------------------------------------------------------------------*
+ * Required types
+ *--------------------------------------------------------------------*/
+
+//--A task for a thread
+
+template <unsigned DIM>
+struct ZoneTask 
+{
+  MZone<DIM> zone;
+  CGNSNameStr zoneName;
+  int zoneIndex;
+  int status;                           // 0 - free
+                                        // 1 - zone processed and waiting in
+                                        //     queue to be written
+                                        // 2 - zone defined and boundaries
+                                        //     processed
+  int indexInOwner;
+  ZoneTask() : status(0), indexInOwner(0) { }
+  void change_status(const int _status) 
+  {
+#pragma omp atomic
+    status = _status;
+  }
+};
+
+//--Information about a zone
+
+struct ZoneInfo 
+{
+  CGNSNameStr name;
+  int cgIndex;
+  ZoneInfo() : cgIndex(-1) { }
+};
+
+template <unsigned DIM>
+int write_CGNS_zones(GModel &model, const int zoneDefinition,
+                     const CGNSOptions &options,
+                     const double scalingFactor, const int vectorDim,
+                     const PhysGroupMap &group,
+                     const int cgIndexFile, const int cgIndexBase)
+{
+
+//--Shared data
+
+  const int numPartitions = model.getNumMeshPartitions();
+  int threadsWorking = omp_get_num_threads();
+                                        // Semaphore for working threads
+  omp_lock_t threadWLock;
+  std::queue<ZoneTask<DIM>*> zoneQueue; // Queue for zones that have been
+                                        // defined and are ready to be written
+  omp_lock_t queueLock;
+  // Next 3 are locked by omp critical
+  int globalZoneIndex = 0;
+  int globalPartition = 1;
+  PhysGroupMap::const_iterator globalItPhysical = group.begin();
+
+//--Initialize omp locks
+
+  omp_init_lock(&threadWLock);
+  omp_init_lock(&queueLock);
+
+//**Spawn threads
+
+  {
+
+//--Master thread (primary task is to define boundary connections and perform
+//--I/O but can also process a zone if idle)
+
+    if(omp_get_thread_num() == 0) {
+      ZoneTask<DIM> zoneTask;
+      MZoneBoundary<DIM> mZoneBoundary;
+      ZoneConnMap zoneConnMap;          // Zone connectivity that is generated
+                                        // whenever a zone is added to
+                                        // mZoneBoundary
+      std::vector<ZoneInfo> zoneInfo; zoneInfo.resize(16);
+      std::vector<double> dBuffer;      // Buffer for double-precision data
+      std::vector<int> iBuffer1, iBuffer2;
+                                        // Buffers for integer data
+      int indexInterface = 0;           // Index for interfaces
+
+      while (threadsWorking || zoneQueue.size()) {
+        if (zoneQueue.size()) {
+
+/*--------------------------------------------------------------------*
+ * Write the zone in the queue
+ *--------------------------------------------------------------------*/
+
+          ZoneTask<DIM> *const writeTask = zoneQueue.front();
+          MZone<DIM> *const writeZone = &writeTask->zone;
+
+//--Write the zone
+
+          // Write the zone node
+          int cgIndexZone;
+          int cgZoneSize[3];
+          cgZoneSize[0] = writeZone->zoneVertVec.size();  // Number of vertices
+          cgZoneSize[1] = writeZone->totalElements();  // Number of elements
+          cgZoneSize[2] = writeZone->numBoVert;  // Number of boundary vertices
+          if(cg_zone_write(cgIndexFile, cgIndexBase,
+                           writeTask->zoneName.c_str(), cgZoneSize,
+                           Unstructured, &cgIndexZone))
+            return cgnsErr();
+//           std::cout << "Zone size: " << cgZoneSize[0] << ' '  // DBG
+//                     << cgZoneSize[1] << ' ' << cgZoneSize[2] << std::endl;
+          // Manually maintain the size of the 'zoneInfo vector'.  'push_back'
+          // is not used because the elements are in order of 'zoneIndex'
+          if(writeTask->zoneIndex >= zoneInfo.size())
+            zoneInfo.resize(std::max(2*static_cast<int>(zoneInfo.size()),
+                                     writeTask->zoneIndex));
+          zoneInfo[writeTask->zoneIndex].name = writeTask->zoneName;
+          zoneInfo[writeTask->zoneIndex].cgIndex = cgIndexZone;
+
+          // Write the grid node
+          int cgIndexGrid;
+          if(cg_grid_write(cgIndexFile, cgIndexBase, cgIndexZone,
+                           "GridCoordinates", &cgIndexGrid))
+            return cgnsErr();
+
+          // Write the grid coordinates
+          int cgIndexCoord;
+          dBuffer.resize(cgZoneSize[0]);
+          // x
+          for(int i = 0; i != cgZoneSize[0]; ++i) 
+              dBuffer[i] = writeZone->zoneVertVec[i]->x()*scalingFactor;
+          if(cg_coord_write(cgIndexFile, cgIndexBase, cgIndexZone, RealDouble,
+                            "CoordinateX", &dBuffer[0], &cgIndexCoord))
+            return cgnsErr();
+          // y
+          for(int i = 0; i != cgZoneSize[0]; ++i)
+            dBuffer[i] = writeZone->zoneVertVec[i]->y()*scalingFactor;
+          if(cg_coord_write(cgIndexFile, cgIndexBase, cgIndexZone, RealDouble,
+                            "CoordinateY", &dBuffer[0], &cgIndexCoord))
+            return cgnsErr();
+          // z
+          if(vectorDim == 3) {
+            for(int i = 0; i != cgZoneSize[0]; ++i)
+              dBuffer[i] = writeZone->zoneVertVec[i]->z()*scalingFactor;
+            if(cg_coord_write(cgIndexFile, cgIndexBase, cgIndexZone, RealDouble,
+                              "CoordinateZ", &dBuffer[0], &cgIndexCoord))
+              return cgnsErr();
+          }
+
+          // Obtain indices sorted in the order the CGNS elements will be
+          // written
+          int indexElemType[MSH_NUM_TYPE];
+          for(int i = 0; i != MSH_NUM_TYPE; ++i) indexElemType[i] = i;
+          std::sort<int*, ElemSortCGNSType>
+            (indexElemType, indexElemType + MSH_NUM_TYPE,
+             ElemSortCGNSType(writeZone->zoneElemConn));
+
+          // Write the element connectivity
+          const int numElemType = writeZone->numElementTypes();
+          int iElemSection = 0;
+          for(int iElemType = 0; iElemType != numElemType; ++iElemType) {
+            const int typeMSHm1 = indexElemType[iElemType];  // typeMSH-1
+            const int typeCGNS = msh2cgns[typeMSHm1][0];
+            const char *elemName;
+            MElement::getInfoMSH(typeMSHm1+1, &elemName);
+            if(typeCGNS == -1) {
+              // This type is not supported in CGNS
+              Message::Warning("Element type %s is not supported in CGNS and "
+                               "has not been written to the file", elemName);
+            }
+            else {
+              int cgIndexSection;
+              //**Replace blanks in 'elemName' with underscore?
+              if(cg_section_write
+                 (cgIndexFile, cgIndexBase, cgIndexZone, elemName,
+                  static_cast<ElementType_t>(typeCGNS), iElemSection + 1,
+                  writeZone->zoneElemConn[typeMSHm1].numElem + iElemSection,
+                  writeZone->zoneElemConn[typeMSHm1].numBoElem + iElemSection,
+                  &writeZone->zoneElemConn[typeMSHm1].connectivity[0],
+                  &cgIndexSection))
+                return cgnsErr();
+              ++iElemSection;
+            }
+          }
+
+//--Process interior connectivity
+
+          mZoneBoundary.interiorBoundaryVertices(writeTask->zoneIndex,
+                                                 *writeZone, zoneConnMap);
+          // Write the connectivity for each zone pair
+          const ZoneConnMap::const_iterator gCEnd = zoneConnMap.end();
+          for(ZoneConnMap::const_iterator gCIt = zoneConnMap.begin();
+              gCIt != gCEnd; ++gCIt) {
+            const int nVert = gCIt->second.vertexPairVec.size();
+            iBuffer1.resize(nVert);
+            iBuffer2.resize(nVert);
+            for(int iVert = 0; iVert != nVert; ++iVert) {
+              const ZoneConnectivity::VertexPair &vp =
+                gCIt->second.vertexPairVec[iVert];
+              iBuffer1[iVert] = vp.vertexIndex1;
+              iBuffer2[iVert] = vp.vertexIndex2;
+            }
+            int cgIndexInterface;
+            // In the first zone
+            if(cg_conn_write
+               (cgIndexFile, cgIndexBase, zoneInfo[gCIt->first.zone1].cgIndex,
+                CGNSNameStr(++indexInterface, "Interface_%d").c_str(), Vertex,
+                Abutting1to1, PointList, nVert, &iBuffer1[0],
+                zoneInfo[gCIt->first.zone2].name.c_str(), Unstructured,
+                PointListDonor, Integer, nVert, &iBuffer2[0],
+                &cgIndexInterface))
+              return cgnsErr();
+            // In the second zone
+            if(cg_conn_write
+               (cgIndexFile, cgIndexBase, zoneInfo[gCIt->first.zone2].cgIndex,
+                CGNSNameStr(++indexInterface, "Interface_%d").c_str(), Vertex,
+                Abutting1to1, PointList, nVert, &iBuffer2[0],
+                zoneInfo[gCIt->first.zone1].name.c_str(), Unstructured,
+                PointListDonor, Integer, nVert, &iBuffer1[0],
+                &cgIndexInterface))
+              return cgnsErr();
           }
-          const int nVPE = (*itElem)->getNumVertices();
-          for(int iVPE = 0; iVPE != nVPE; ++iVPE)
-            elementConn[iElemData++] = (*itElem)->getVertex(iVPE)->getNum();
+
+//--Pop from the queue
+
+          omp_set_lock(&queueLock);
+          zoneQueue.pop();
+          omp_unset_lock(&queueLock);
+
+//--Task finished
+
+          writeTask->change_status(0);
         }
-      }
-      else {
-        const int nVPE = numVPE;
-        for(InteriorElemSet::const_iterator itElem = interiorElem.begin();
-            itElem != itEndIElem; ++itElem)
-          for(int iVPE = 0; iVPE != nVPE; ++iVPE)
-            elementConn[iElemData++] = (*itElem)->getVertex(iVPE)->getNum();
-      }
 
-      // Write to CGNS file
-      int cgIndexSection;
-      if(cg_section_write(cgIndexFile, cgIndexBase, cgIndexZone, "Section_Main",
-                          cgElemType, 1, cgZoneSize[1], boundaryElem.size(),
-                          &elementConn[0], &cgIndexSection)) return cgnsErr();
-    }
+/*--------------------------------------------------------------------*
+ * No zones waiting in the queue, process a zone
+ *--------------------------------------------------------------------*/
 
-    meshNeighbours.clear();
-    ++indexZone;
+        else {
+          PhysGroupMap::const_iterator itPhysicalBegin;
+          PhysGroupMap::const_iterator itPhysicalEnd;
+          int zoneIndex;
+          int partition;
+          if(get_zone_definition
+             (model, zoneDefinition, options, numPartitions, group,
+              globalZoneIndex, globalPartition, globalItPhysical,
+              zoneTask.zoneIndex, partition, itPhysicalBegin, itPhysicalEnd,
+              zoneTask.zoneName)) {
+            omp_set_lock(&threadWLock);
+            --threadsWorking;
+            omp_unset_lock(&threadWLock);
+          }
+          else {
+            zoneTask.zone.clear();
+            // Loop through all physical in the zone definition
+            for(PhysGroupMap::const_iterator itPhysical = itPhysicalBegin;
+                itPhysical != itPhysicalEnd; ++itPhysical) {
+              // Add elements from all entities in the physical with defined
+              // partition number
+              zoneTask.zone.template add_elements_in_entities
+                <typename DimTr<DIM>::EntityT,
+                typename std::vector<GEntity*>::const_iterator>
+                (itPhysical->second.begin(), itPhysical->second.end(),
+                 partition);
+            }
+            // Process the zone
+            zoneTask.zone.zoneData();
+            translateElementMSH2CGNS(zoneTask.zone.zoneElemConn);
+            // Add to the queue to get written
+            zoneTask.change_status(1);
+            omp_set_lock(&queueLock);
+            zoneQueue.push(&zoneTask);
+            omp_unset_lock(&queueLock);
+          }
+        }
+      }  // End master thread loop
 
-  }
+/*--------------------------------------------------------------------*
+ * Write the remaining unconnected vertices as boundary conditions
+ *--------------------------------------------------------------------*/
 
-//--Close the file
+      ZoneBoVec zoneBoVec;              // from 'MZoneBoundary.h'
+      if(mZoneBoundary.exteriorBoundaryVertices(zoneBoVec) == 0) {  // Have BC
+
+        // Sort by zone index and then entity index
+        const int numBoVert = zoneBoVec.size();
+        std::vector<int> iZBV(numBoVert);
+        for(int i = 0; i != numBoVert; ++i) iZBV[i] = i;
+        std::sort<int*, ZoneBoVecSort>(&iZBV[0], &iZBV[numBoVert],
+                                       ZoneBoVecSort(zoneBoVec));
+
+        dBuffer.reserve(1024);
+        iBuffer1.reserve(1024);
+
+        int iVert = 0;
+        while(iVert != numBoVert) {
+          dBuffer.clear();
+          iBuffer1.clear();
+          const int zoneIndex = zoneBoVec[iZBV[iVert]].zoneIndex;
+          const int patchIndex = zoneBoVec[iZBV[iVert]].bcPatchIndex;
+          const int iVertStart = iVert;
+          while(iVert != numBoVert &&
+                zoneBoVec[iZBV[iVert]].zoneIndex == zoneIndex &&
+                zoneBoVec[iZBV[iVert]].bcPatchIndex == patchIndex) {
+            VertexBoundary &vertBo = zoneBoVec[iZBV[iVert]];
+            dBuffer.push_back(vertBo.normal[0]);
+            dBuffer.push_back(vertBo.normal[1]);
+            if(vectorDim == 3) dBuffer.push_back(vertBo.normal[2]);
+            iBuffer1.push_back(vertBo.vertexIndex);
+            ++iVert;
+          }
+          const int numBCVert = iVert - iVertStart;
+
+          int cgIndexBoco;
+          if(cg_boco_write(cgIndexFile, cgIndexBase,
+                           zoneInfo[zoneIndex].cgIndex,
+                           CGNSNameStr(patchIndex, "Patch %d").c_str(),
+                           BCTypeNull, PointList, numBCVert, &iBuffer1[0],
+                           &cgIndexBoco))
+            return cgnsErr();
+
+          int normalIndex;
+          if(cg_boco_normal_write(cgIndexFile, cgIndexBase,
+                                  zoneInfo[zoneIndex].cgIndex, cgIndexBoco,
+                                  &normalIndex, 1, RealDouble, &dBuffer[0]))
+            return cgnsErr();
+        }
+      }
 
-  if(cg_close(cgIndexFile)) return cgnsErr();
+//       std::cout << "Leaving master thread\n";  // DBG
+    }  // End master thread instructions
+
+  }  // End omp parallel section
+
+//--Destroy omp locks
+
+  omp_destroy_lock(&threadWLock);
+  omp_destroy_lock(&queueLock);
 
-  std::cout << "Done\n";
-  return 1;
 }
 
 #else
 
 int GModel::readCGNS(const std::string &name)
 {
-  Msg::Error("This version of Gmsh was compiled without CGNS support");
+  Message::Error("This version of Gmsh was compiled without CGNS support");
   return 0;
 }
 
-int GModel::writeCGNS(const std::string &name, double scalingFactor)
+int GModel::writeCGNS(const std::string &name, const int zoneDefinition,
+                      const CGNSOptions &options, double scalingFactor)
 {
-  Msg::Error("This version of Gmsh was compiled without CGNS support");
+  Message::Error("This version of Gmsh was compiled without CGNS support");
   return 0;
 }
 
 #endif
-
diff --git a/Geo/MEdge.h b/Geo/MEdge.h
index 351e4c1456..d93d5297a8 100644
--- a/Geo/MEdge.h
+++ b/Geo/MEdge.h
@@ -34,6 +34,12 @@ class MEdge {
   inline MVertex *getSortedVertex(const int i) const { return _v[int(_si[i])]; }
   inline MVertex *getMinVertex() const { return _v[int(_si[0])]; }
   inline MVertex *getMaxVertex() const { return _v[int(_si[1])]; }
+  SVector3 scaledTangent() const
+  {
+    return SVector3(_v[1]->x() - _v[0]->x(), 
+                    _v[1]->y() - _v[0]->y(),
+                    _v[1]->z() - _v[0]->z());
+  }
   SVector3 tangent() const
   {
     SVector3 t(_v[1]->x() - _v[0]->x(), 
@@ -68,6 +74,12 @@ class MEdge {
   }
 };
 
+inline bool operator==(const MEdge &e1, const MEdge &e2)
+{
+  return (e1.getMinVertex() == e2.getMinVertex() &&
+          e1.getMaxVertex() == e2.getMaxVertex());
+}
+
 inline bool operator!=(const MEdge &e1, const MEdge &e2)
 {
   return (e1.getMinVertex() != e2.getMinVertex() ||
@@ -77,8 +89,7 @@ inline bool operator!=(const MEdge &e1, const MEdge &e2)
 struct Equal_Edge : public std::binary_function<MEdge, MEdge, bool> {
   bool operator()(const MEdge &e1, const MEdge &e2) const
   {
-    return (e1.getMinVertex() == e2.getMinVertex() &&
-            e1.getMaxVertex() == e2.getMaxVertex());
+    return (e1 == e2);
   }
 };
 
diff --git a/Geo/MFace.h b/Geo/MFace.h
index d00b22d598..9e7eb86730 100644
--- a/Geo/MFace.h
+++ b/Geo/MFace.h
@@ -91,13 +91,26 @@ class MFace {
   }
 };
 
+inline bool operator==(const MFace &f1, const MFace &f2)
+{
+  return (f1.getSortedVertex(0) == f2.getSortedVertex(0) &&
+          f1.getSortedVertex(1) == f2.getSortedVertex(1) &&
+          f1.getSortedVertex(2) == f2.getSortedVertex(2) &&
+          f1.getSortedVertex(3) == f2.getSortedVertex(3));
+}
+
+inline bool operator!=(const MFace &f1, const MFace &f2)
+{
+  return (f1.getSortedVertex(0) != f2.getSortedVertex(0) ||
+          f1.getSortedVertex(1) != f2.getSortedVertex(1) ||
+          f1.getSortedVertex(2) != f2.getSortedVertex(2) ||
+          f1.getSortedVertex(3) != f2.getSortedVertex(3));
+}
+
 struct Equal_Face : public std::binary_function<MFace, MFace, bool> {
   bool operator()(const MFace &f1, const MFace &f2) const
   {
-    return (f1.getSortedVertex(0) == f2.getSortedVertex(0) &&
-            f1.getSortedVertex(1) == f2.getSortedVertex(1) &&
-            f1.getSortedVertex(2) == f2.getSortedVertex(2) &&
-            f1.getSortedVertex(3) == f2.getSortedVertex(3));
+    return (f1 == f2);
   }
 };
 
diff --git a/Geo/MNeighbour.h b/Geo/MNeighbour.h
deleted file mode 100644
index c2bb20f8ab..0000000000
--- a/Geo/MNeighbour.h
+++ /dev/null
@@ -1,1941 +0,0 @@
-#ifndef _MNEIGHBOUR_H_
-#define _MNEIGHBOUR_H_
-
-// Copyright (C) 2006 S. Guzik, C. Geuzaine, J.-F. Remacle
-//
-// This program is free software; you can redistribute it and/or modify
-// it under the terms of the GNU General Public License as published by
-// the Free Software Foundation; either version 2 of the License, or
-// (at your option) any later version.
-//
-// This program is distributed in the hope that it will be useful,
-// but WITHOUT ANY WARRANTY; without even the implied warranty of
-// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-// GNU General Public License for more details.
-//
-// You should have received a copy of the GNU General Public License
-// along with this program; if not, write to the Free Software
-// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
-// USA.
-// 
-// Please report all bugs and problems to <gmsh@geuz.org>.
-
-/*******************************************************************************
- *
- * - The classes in this file construct a database of the elements that share
- *   the lower-dimensional bounding objects (e.g., vertex, edge, or face).
- *   These lower-dimensional objects are referred to in general as polytopes.
- * - Templated traits classes are used extensively to define the characteristics
- *   of the elements and the entities.
- *
- ******************************************************************************/
-
-#include <algorithm>
-#include <iterator>
-#include <vector>
-#include "MElement.h"
-#include "GEdge.h"
-#include "GFace.h"
-#include "GRegion.h"
-#include "MEdgeHash.h"
-#include "MFaceHash.h"
-#include "GmshDefines.h"
-
-// #define HAVE_HASH_MAP
-
-#if defined(HAVE_HASH_MAP)
-#include "HashMap.h"
-#endif
-
-
-/*==============================================================================
- * File scope types
- *============================================================================*/
-
-namespace {
-
-typedef std::list<MElement*> Neighbours;
-typedef Neighbours::const_iterator NeighboursConstIterator;
-typedef Neighbours::iterator NeighboursIterator;
-
-struct Range_t
-{
-  int num;
-  NeighboursIterator begin;
-  Range_t() : num(0) { }
-};
-
-//--Use a hash map for neighbour lookup if possible, otherwise a map will do
-
-#if defined(HAVE_HASH_MAP)
-struct Hash_VertexPtr : public std::unary_function<MVertex*, size_t>
-{
-  size_t operator()(const MVertex *const v) const
-  {
-    return HashFNV1a<sizeof(MVertex*)>::eval(v);
-  }
-};
-typedef HASH_MAP<MVertex*, Range_t, Hash_VertexPtr, std::equal_to<MVertex*> >
-VertNRange;
-typedef HASH_MAP<MEdge, Range_t, Hash_Edge, Equal_Edge> EdgeNRange;
-typedef HASH_MAP<MFace, Range_t, Hash_Face, Equal_Face> FaceNRange;
-#else
-typedef std::map<MVertex*, Range_t, std::less<MVertex*> > VertNRange;
-typedef std::map<MEdge, Range_t, Less_Edge> EdgeNRange;
-typedef std::map<MFace, Range_t, Less_Face> FaceNRange;
-#endif
-typedef VertNRange::iterator VertNRangeIterator;
-typedef VertNRange::const_iterator VertNRangeConstIterator;
-typedef EdgeNRange::iterator EdgeNRangeIterator;
-typedef EdgeNRange::const_iterator EdgeNRangeConstIterator;
-typedef FaceNRange::iterator FaceNRangeIterator;
-typedef FaceNRange::const_iterator FaceNRangeConstIterator;
-
-
-/*==============================================================================
- * Traits classes - that just perform compile-time checks for valid argument
- * types
- *============================================================================*/
-
-//--This traits class checks for a pointer and strips it
-
-// NOTE:  If the compiler sent you here, you should be using a pointer.
-template <typename T> struct StripPointer;
-template <typename T> struct StripPointer<T*> { typedef T Type; };
-    
-//--This traits class checks for valid types of entities and iterators when the
-//--entity is explicitly specified
-
-// NOTE:  If the compiler sent you here, you're using an invalid entity and/or
-// invalid iterator.  Valid pairs are shown below
-template <typename Ent, typename Iter> struct ValidSpecEntityIterator;
-template <> struct ValidSpecEntityIterator<GEdge, GEdge*>
-{ typedef void Type; };
-template <> struct ValidSpecEntityIterator<GEdge, GEntity*>
-{ typedef void Type; };
-template <> struct ValidSpecEntityIterator<GFace, GFace*>
-{ typedef void Type; };
-template <> struct ValidSpecEntityIterator<GFace, GEntity*>
-{ typedef void Type; };
-template <> struct ValidSpecEntityIterator<GRegion, GRegion*>
-{ typedef void Type; };
-template <> struct ValidSpecEntityIterator<GRegion, GEntity*>
-{ typedef void Type; };
-
-//--This traits class checks for valid types of entities and iterators when the
-//--entity is not specified.  It also strips the pointer.
-
-// NOTE:  If the compiler sent you here, you're using an invalid entity or
-// iterator of entities
-template <typename Iter> struct ValidEntityIterator;
-template <> struct ValidEntityIterator<GEdge*>
-{ typedef GEdge Type; };
-template <> struct ValidEntityIterator<GFace*>
-{ typedef GFace Type; };
-template <> struct ValidEntityIterator<GRegion*>
-{ typedef GRegion Type; };
-
-//--This traits class checks for valid element types
-
-// NOTE:  If the compiler sent you here, you're using an invalid iterator for an
-// element.  What's listed here are NOT the only valid types.  Any iterator with
-// a value type of a pointer to an element should work.
-template <typename Elem_o1> struct ValidElement;
-template <> struct ValidElement<MElement> { typedef void Type; };
-template <> struct ValidElement<MLine> { typedef void Type; };
-template <> struct ValidElement<MTriangle> { typedef void Type; };
-template <> struct ValidElement<MQuadrangle> { typedef void Type; };
-template <> struct ValidElement<MTetrahedron> { typedef void Type; };
-template <> struct ValidElement<MHexahedron> { typedef void Type; };
-template <> struct ValidElement<MPrism> { typedef void Type; };
-template <> struct ValidElement<MPyramid> { typedef void Type; };
-
-
-/*==============================================================================
- * Traits classes - that return information about a type
- *============================================================================*/
-
-//--This traits class gives the first-order base type for a second-order element
-
-template <typename Order2> struct ElemBaseOrder
-{ typedef Order2 Order1; };
-template <> struct ElemBaseOrder<MLine3>
-{ typedef MLine Order1; };
-template <> struct ElemBaseOrder<MTriangle6>
-{ typedef MTriangle Order1; };
-template <> struct ElemBaseOrder<MQuadrangle8>
-{ typedef MQuadrangle Order1; };
-template <> struct ElemBaseOrder<MQuadrangle9>
-{ typedef MQuadrangle Order1; };
-template <> struct ElemBaseOrder<MTetrahedron10>
-{ typedef MTetrahedron Order1; };
-template <> struct ElemBaseOrder<MHexahedron20>
-{ typedef MHexahedron Order1; };
-template <> struct ElemBaseOrder<MHexahedron27>
-{ typedef MHexahedron Order1; };
-template <> struct ElemBaseOrder<MPrism15>
-{ typedef MPrism Order1; };
-template <> struct ElemBaseOrder<MPrism18>
-{ typedef MPrism Order1; };
-template <> struct ElemBaseOrder<MPyramid13>
-{ typedef MPyramid Order1; };
-template <> struct ElemBaseOrder<MPyramid14>
-{ typedef MPyramid Order1; };
-
-//--This traits class describes describes the lower dimensional bounding
-//--polytopes of a first-order element
-
-template <typename Elem> struct ElemTr;
-template <> struct ElemTr<MElement>
-{ enum { numVertex = 0, numEdge = 0, numFace = 0 }; };
-template <> struct ElemTr<MLine>
-{ enum { numVertex = 2, numEdge = 0, numFace = 0 }; };
-template <> struct ElemTr<MTriangle>
-{ enum { numVertex = 3, numEdge = 3, numFace = 0 }; };
-template <> struct ElemTr<MQuadrangle>
-{ enum { numVertex = 4, numEdge = 4, numFace = 0 }; };
-template <> struct ElemTr<MTetrahedron>
-{ enum { numVertex = 4, numEdge = 6, numFace = 4 }; };
-template <> struct ElemTr<MHexahedron>
-{ enum { numVertex = 8, numEdge = 12, numFace = 6 }; };
-template <> struct ElemTr<MPrism>
-{ enum { numVertex = 6, numEdge = 9, numFace = 5 }; };
-template <> struct ElemTr<MPyramid>
-{ enum { numVertex = 5, numEdge = 8, numFace = 5 }; };
-
-//--This traits class gives the number of element types in entity Ent
-
-template <typename Ent> struct EntTr;
-template <> struct EntTr<GEdge> { enum { numElemTypes = 1 }; };
-template <> struct EntTr<GFace> { enum { numElemTypes = 2 }; };
-template <> struct EntTr<GRegion> { enum { numElemTypes = 4 }; };
-
-//--This traits class gives iterator types and begin and end iterators for
-//--element type number N in entity Ent.
-
-template <typename Ent, int N> struct EntElemTr;
-template <> struct EntElemTr<GEdge, 1> {
-  typedef MLine Elem;
-  typedef std::vector<MLine*>::const_iterator ConstElementIterator;
-  typedef std::vector<MLine*>::iterator ElementIterator;
-  static ConstElementIterator begin(const GEdge *const edge)
-  {
-    return edge->lines.begin();
-  }
-  static ConstElementIterator end(const GEdge *const edge)
-  {
-    return edge->lines.end();
-  }
-};
-template <> struct EntElemTr<GFace, 1> {
-  typedef MQuadrangle Elem;
-  typedef std::vector<MQuadrangle*>::const_iterator ConstElementIterator;
-  typedef std::vector<MQuadrangle*>::iterator ElementIterator;
-  static ConstElementIterator begin(const GFace *const face)
-  {
-    return face->quadrangles.begin();
-  }
-  static ConstElementIterator end(const GFace *const face)
-  {
-    return face->quadrangles.end();
-  }
-};
-template <> struct EntElemTr<GFace, 2> {
-  typedef MTriangle Elem;
-  typedef std::vector<MTriangle*>::const_iterator ConstElementIterator;
-  typedef std::vector<MTriangle*>::iterator ElementIterator;
-  static ConstElementIterator begin(const GFace *const face)
-  {
-    return face->triangles.begin();
-  }
-  static ConstElementIterator end(const GFace *const face)
-  {
-    return face->triangles.end();
-  }
-};
-template <> struct EntElemTr<GRegion, 1> {
-  typedef MPyramid Elem;
-  typedef std::vector<MPyramid*>::const_iterator ConstElementIterator;
-  typedef std::vector<MPyramid*>::iterator ElementIterator;
-  static ConstElementIterator begin(const GRegion *const region)
-  {
-    return region->pyramids.begin();
-  }
-  static ConstElementIterator end(const GRegion *const region)
-  {
-    return region->pyramids.end();
-  }
-};
-template <> struct EntElemTr<GRegion, 2> {
-  typedef MPrism Elem;
-  typedef std::vector<MPrism*>::const_iterator ConstElementIterator;
-  typedef std::vector<MPrism*>::iterator ElementIterator;
-  static ConstElementIterator begin(const GRegion *const region)
-  {
-    return region->prisms.begin();
-  }
-  static ConstElementIterator end(const GRegion *const region)
-  {
-    return region->prisms.end();
-  }
-};
-template <> struct EntElemTr<GRegion, 3> {
-  typedef MHexahedron Elem;
-  typedef std::vector<MHexahedron*>::const_iterator ConstElementIterator;
-  typedef std::vector<MHexahedron*>::iterator ElementIterator;
-  static ConstElementIterator begin(const GRegion *const region)
-  {
-    return region->hexahedra.begin();
-  }
-  static ConstElementIterator end(const GRegion *const region)
-  {
-    return region->hexahedra.end();
-  }
-};
-template <> struct EntElemTr<GRegion, 4> {
-  typedef MTetrahedron Elem;
-  typedef std::vector<MTetrahedron*>::const_iterator ConstElementIterator;
-  typedef std::vector<MTetrahedron*>::iterator ElementIterator;
-  static ConstElementIterator begin(const GRegion *const region)
-  {
-    return region->tetrahedra.begin();
-  }
-  static ConstElementIterator end(const GRegion *const region)
-  {
-    return region->tetrahedra.end();
-  }
-};
-
-//--This is a traits/policy class for the lower-dimension polytopes that bound
-//--an element (i.e., vertex, edge, or face).  It returns the corresponding
-//--range type and extracts from the element the function to get the polytope.
-//--Note that MVertex is a pointer because it exists somewhere whereas MEdge and
-//--MFace are constructed using MVertex's by getEdge() and getFace().
-
-template <typename Polytope> struct PolytopeTr;
-template <> struct PolytopeTr<MVertex*> {
-  typedef VertNRange PolytopeNRange;
-  typedef VertNRangeConstIterator PNRConstIterator;
-  template <typename Elem>
-  static MVertex *getPolytope(Elem *const element, const int nPolytope)
-  {
-    return element->getVertex(nPolytope);
-  }
-};
-template <> struct PolytopeTr<MEdge> {
-  typedef EdgeNRange PolytopeNRange;
-  typedef EdgeNRangeConstIterator PNRConstIterator;
-  template <typename Elem>
-  static MEdge getPolytope(Elem *const element, const int nPolytope)
-  {
-    return element->getEdge(nPolytope);
-  }
-};
-template <> struct PolytopeTr<MFace> {
-  typedef FaceNRange PolytopeNRange;
-  typedef FaceNRangeConstIterator PNRConstIterator;
-  template <typename Elem>
-  static MFace getPolytope(Elem *const element, const int nPolytope)
-  {
-    return element->getFace(nPolytope);
-  }
-};
-
-}  // End of unnamed namespace
-
-
-/*******************************************************************************
- *
- * class: PolytopeIterator
- *
- * Purpose
- * =======
- *
- *   Provides an iterator to iterate through the unique vertices, edges and,
- *   faces defined by the (hash_)maps.  An iterator to the (hash_)map has a
- *   value type of a (key, value) pair.  This iterator makes the container look
- *   like it only contains the bounding polytope.
- *
- * Constructors
- * ============
- *
- *   The class takes one template argument <Polytope> which should be either
- *   MVertex*, MEdge, or MFace.
- *
- *   PolytopeIterator()
- *                     -- default constructor
- *
- *   PolytopeIterator(const BaseIterator &baseIter_in)
- *                     -- the base iterator is the real iterator to the
- *                        (hash_)map
- *
- *   PolytopeIterator(const  PolytopeIterator &iter)
- *                     -- copy
- *
- * Destructor
- * ==========
- *
- *   ~PolytopeIterator -- synthesized
- *
- * Member Functions
- * ================
- *
- *   int num_neighbours()
- *                     -- returns number of elements sharing the polytope
- *
- * Operators
- * =========
- *
- *   Includes typical bidirectional iterator operators {=, ==, !=, *, ->, ++(),
- *   --(), ()++, ()--} with * and -> dereferencing to the polytope.
- *
- * Notes
- * =====
- *
- *   - Only constant iterators are supported.
- *   
- ******************************************************************************/
-
-template<typename Polytope>
-class PolytopeIterator  // Actually only a const_iterator
-{
-
- public:
-
-  // The base iterator iterates through the (hash_)maps defining the range of
-  // element neighbours for a polytope
-  typedef typename PolytopeTr<Polytope>::PNRConstIterator BaseIterator;
-  typedef PolytopeIterator Self;
-
-  // Standard traits
-  typedef std::bidirectional_iterator_tag iterator_category;
-  typedef Polytope value_type;
-  typedef BaseIterator difference_type;
-  typedef const Polytope* pointer;
-  typedef const Polytope& reference;
-
-//--Constructors
-
-  PolytopeIterator() : baseIter() { }
-  PolytopeIterator(const BaseIterator &baseIter_in) : baseIter(baseIter_in) { }
-  PolytopeIterator(const PolytopeIterator &iter) : baseIter(iter.baseIter) { }
-  PolytopeIterator& operator=(const PolytopeIterator& iter)
-  {
-    if(&iter != this) baseIter = iter.baseIter;
-    return *this;
-  }
-
-//--Comparison
-
-  bool operator==(const Self& iter) const { return baseIter == iter.baseIter; }
-  bool operator!=(const Self& iter) const { return baseIter != iter.baseIter; }
-
-//--Dereference and arrow operator
-
-  reference operator*() const { return baseIter->first; }
-  pointer operator->() const { return &baseIter->first; }
-
-//--Increment
-
-  // Prefix
-  Self& operator++()
-  {
-    ++baseIter;
-    return *this;
-  }
-
-  // Postfix
-  Self operator++(int)
-  {
-    Self tmp = *this;
-    ++*this;
-    return tmp;
-  }
-
-//--Decrement
-
-  // Prefix
-  Self& operator--()
-  {
-    --baseIter;
-    return *this;
-  }
-
-  // Postfix
-  Self operator--(int)
-  {
-    Self tmp = *this;
-    --*this;
-    return tmp;
-  }
-
-//--Special
-
-  int num_neighbours() { return baseIter->second.num; }
-
- private:
-  NeighboursConstIterator begin_neighbours() { return baseIter->second.begin; }
-
-//--Date members
-
-  BaseIterator baseIter;
-
-//--Friends
-
-  friend class MNeighbour;
-
-};
-
-
-/*******************************************************************************
- *
- * class: MNeighbour
- *
- * Purpose
- * =======
- *
- *   Generates a database for quick lookup of elements sharing a
- *   MVertex*, MEdge, or MFace (referred to in general as polytopes.  I.e., a
- *   geometrical construct bounding the element with a lower dimension than the
- *   element).  The neighbour elements are stored in a list.  Maps or hash maps
- *   provide an index into the list for a given polytope.  In addition to the
- *   element neighbours, this class provides iterators to the unique vertices,
- *   edges, and faces in the database.
- *
- * Iterators
- * =========
- *
- *   Vertex_const_iterator
- *                     -- iterator to sequence of unique vertices
- *   Edge_const_iterator
- *                     -- iterator to sequence of unique edges
- *   Face_const_iterator
- *                     -- iterator to sequence of unique faces
- *
- * Constructors
- * ============
- *
- *   MNeighbour()      -- default constructor does nothing
- *
- * Destructor
- * ==========
- *
- *   ~MNeighbour()     -- synthesized
- *
- * Member Functions
- * ================
- *
- * Iterators
- * ---------
- *
- *   Vertex_const_iterator vertex_begin()
- *                     -- first unique vertex
- *
- *   Vertex_const_iterator vertex_end()
- *                     -- one-past-last unique vertex
- *
- *   Edeg_const_iterator edge_begin()
- *                     -- first unique edge
- *
- *   Edge_const_iterator edge_end()
- *                     -- one-past-last unique edge
- *
- *   Face_const_iterator face_begin()
- *                     -- first unique face
- *
- *   Face_const_iterator face_end()
- *                     -- one-past-last unique face
- *
- * Routines for adding elements
- * ----------------------------
- *
- *   The routines that generate a database of neighbours from container of
- *   entities (add_elements_in_entitie) require knowledge of the specific type
- *   of entity.  In cases where a container contains specific entities, the type
- *   need not be stated.  In cases where the container contains GEntity*, the
- *   user must state the actual entity type.  Extensive compile-time type
- *   checking is employed to make sure the routines are used correctly.
- *
- *   template <typename Ent, typename EntIter>
- *   void add_elements_in_entities(EntIter begin, EntIter end)
- *                     -- adds elements in entities to neighbours database.
- *                        Entity type <Ent> must be explicitly specified.  The
- *                        value type of the iterators must be a pointer to an
- *                        entity.  It is expected, but not required to be a
- *                        pointer to a GEntity.
- *     begin              (I) iterator/pointer to first entity pointer
- *     end                (I) iterator/pointer to one-past-last entity pointer
- *
- *   template <typename Ent, typename EntP>
- *   void add_elements_in_entities(EntP **pEnt, unsigned int n)
- *                     -- adds elements in entities to neighbours database.
- *                        Entity type <Ent> must be explicitly specified.  The
- *                        array element must be a pointer to an element.  It is
- *                        expected, but not required to be a pointer to a
- *                        GEntity.
- *     pEnt               (I) first entity in the array
- *     n                  (I) number of entities in the array
- *
- *   template <typename EntIter>
- *   void add_elements_in_entities(EntIter begin, EntIter end)
- *                     -- adds elements in entities to neighbours database.
- *                        The entity type must not be GEntity.  The value type
- *                        of the iterators must be a pointer to an entity.
- *     begin              (I) iterator/pointer to first entity pointer
- *     end                (I) iterator/pointer to one-past-last entity pointer
- *
- *   template <typename Ent>
- *   void add_elements_in_entities(Ent **pEnt, unsigned int n)
- *                     -- adds elements in entities to neighbours database.
- *                        The entity type must not be GEntity.  The array
- *                        element must be a pointer to an element.
- *     pEnt               (I) first entity in the array
- *     n                  (I) number of entities in the array
- *
- *   template <typename EntPtr>
- *   void add_elements_in_entity(EntPtr entity)
- *                     -- adds elements in a single entity to the neighbours
- *                        database.  The entity type must not be GEntity.  The
- *                        dereferenced quantity must be a pointer to an entity.
- *     entity             (I) pointer to an entity
- *
- *   The routines that generate a database of neighbours from containers of
- *   elements do not require the user to state the specific type of element.
- *   If the container constains a specific element pointer, compiler-time
- *   polymorphism is used.  If the container contains the base pointers
- *   MElement*, run-time switching based on element introspection is used.
- *
- *   template <typename ElemIter>
- *   void add_elements(ElemIter begin, ElemIter end)
- *                     -- adds elements from iterator sequence to the neighbours
- *                        database.
- *     begin              (I) iterator/pointer to first element pointer
- *     end                (I) iterator/pointer to one-past-last element pointer
- *
- *   template <typename Elem>
- *   void add_elements(Elem **pElem, const unsigned int n)
- *                     -- adds elements in an array to the neighbours database.
- *     pElem              (I) first element in the array
- *     n                  (I) number of elements in the array
- *
- *   void clear()      -- clears the neighbours database.  Use this before
- *                        constructing a new database of element neighbours.
- *                        Otherwise, new elements are added to the existing
- *                        database.
- *
- * Routines for querying the neighbours database
- * ---------------------------------------------
- *
- *   int max_vertex_neighbours()
- *                     -- returns max number of elements sharing a vertex.  This
- *                        is also the max number of elements sharing any
- *                        bounding polytope.
- *
- *   int max_edge_neighbours()
- *                     -- returns max number of elements sharing an edge.
- *
- *   int max_face_neighbours()
- *                     -- returns max number of elements sharing a face.
- *
- *   int vertex_neighbours(Vertex_const_iterator &itVert,
- *                          std::vector<MElement*> &vElem)
- *                     -- return all elements sharing a vertex
- *     itVert             (I) iterator to a unique vertex.  Avoids a find()
- *     vElem              (O) vector with elements loaded via push_back
- *     return             number of neighbour elements
- *
- *   int vertex_neighbours(Vertex_const_iterator &itVert, MElement **pElem)
- *                     -- return all elements sharing a vertex
- *     itVert             (I) iterator to a unique vertex.  Avoids a find()
- *     pElem              (O) array with elements loaded starting at pElem[0]
- *     return             number of neighbour elements
- *
- *   int vertex_neighbours(MVertex *const vertex, std::vector<MElement*> &vElem)
- *                     -- returns all elements sharing a vertex
- *     vertex             (I) pointer to vertex.  It must be found in (hash_)map
- *     vElem              (O) vector with elements loaded via push_back
- *     return             number of neighbour elements.  0 if vertex not found
- *
- *   int vertex_neighbours(MVertex *const vertex, MElement **pElem)
- *                     -- returns all elements sharing a vertex
- *     vertex             (I) pointer to vertex.  It must be found in (hash_)map
- *     pElem              (O) array with elements loaded starting at pElem[0]
- *     return             number of neighbour elements.  0 if vertex not found
- *
- *   int edge_neighbours(Edge_const_iterator &itEdge,
- *                        std::vector<MElement*> &vElem)
- *                     -- return all elements sharing a edge
- *     itEdge             (I) iterator to a unique edge.  Avoids a find()
- *     vElem              (O) vector with elements loaded via push_back
- *     return             number of neighbour elements
- *
- *   int edge_neighbours(Edge_const_iterator &itEdge, MElement **pElem)
- *                     -- return all elements sharing a edge
- *     itEdge             (I) iterator to a unique edge.  Avoids a find()
- *     pElem              (O) array with elements loaded starting at pElem[0]
- *     return             number of neighbour elements
- *
- *   int edge_neighbours(MEdge *const edge, std::vector<MElement*> &vElem)
- *                     -- returns all elements sharing a edge
- *     edge               (I) pointer to edge.  It must be found in (hash_)map
- *     vElem              (O) vector with elements loaded via push_back
- *     return             number of neighbour elements.  0 if edge not found
- *
- *   int edge_neighbours(MEdge *const edge, MElement **pElem)
- *                     -- returns all elements sharing a edge
- *     edge               (I) pointer to edge.  It must be found in (hash_)map
- *     pElem              (O) array with elements loaded starting at pElem[0]
- *     return             number of neighbour elements.  0 if edge not found
- *
- *   int face_neighbours(Face_const_iterator &itFace,
- *                        std::vector<MElement*> &vElem)
- *                     -- return all elements sharing a face
- *     itFace             (I) iterator to a unique face.  Avoids a find()
- *     vElem              (O) vector with elements loaded via push_back
- *     return             number of neighbour elements
- *
- *   int face_neighbours(Face_const_iterator &itFace, MElement **pElem)
- *                     -- return all elements sharing a face
- *     itFace             (I) iterator to a unique face.  Avoids a find()
- *     pElem              (O) array with elements loaded starting at pElem[0]
- *     return             number of neighbour elements
- *
- *   int face_neighbours(MFace *const face, std::vector<MElement*> &vElem)
- *                     -- returns all elements sharing a face
- *     face               (I) pointer to face.  It must be found in (hash_)map
- *     vElem              (O) vector with elements loaded via push_back
- *     return             number of neighbour elements.  0 if face not found
- *
- *   int face_neighbours(MFace *const face, MElement **pElem)
- *                     -- returns all elements sharing a face
- *     face               (I) pointer to face.  It must be found in (hash_)map
- *     pElem              (O) array with elements loaded starting at pElem[0]
- *     return             number of neighbour elements.  0 if face not found
- *
- *   int vertex_num_neighbours(MVertex *const vertex)
- *                     -- returns the number of elements sharing a vertex
- *
- *   int edge_num_neighbours(const MEdge& edge)
- *                     -- returns the number of elements sharing an edge
- *
- *   int face_num_neighbours(const MFace& face)
- *                     -- returns the number of elements sharing a face
- *
- *   template <typename Elem>
- *   int all_element_neighbours(Elem *const element,
- *                              std::vector<MElement*> &vElem,
- *                              const bool exclusive = true)
- *                     -- find all the elements connected to the given element
- *     element            (I) element to find neighbours to
- *     vElem              (O) vector with elements loaded via push_back
- *     exclusive          (I) T - exclude the given element in vElem
- *                            F - include the given element in vElem
- *     return             number of elements added to vElem
- *
- *   template <typename Elem>
- *   int all_element_neighbours(Elem *const element, MElement **pElem
- *                              const bool exclusive = true)
- *                     -- find all the elements connected to the given element
- *     element            (I) element to find neighbours to
- *     pElem              (O) array with elements loaded starting at pElem[0]
- *     exclusive          (I) T - exclude the given element in vElem
- *                            F - include the given element in vElem
- *     return             number of elements added to pElem
- *
- *   template <typename Elem>
- *   int all_element_edge_neighbours(Elem *const element,
- *                                   std::vector<MElement*> &vElem,
- *                                   const bool exclusive = true)
- *                     -- find all the elements connected to the edges of a
- *                        given element
- *     element            (I) element to find neighbours to
- *     vElem              (O) vector with elements loaded via push_back
- *     exclusive          (I) T - exclude the given element in vElem
- *                            F - include the given element in vElem
- *     return             number of elements added to vElem
- *
- *   template <typename Elem>
- *   int all_element_edge_neighbours(Elem *const element, MElement **pElem
- *                                   const bool exclusive = true)
- *                     -- find all the elements connected to the edges of a
- *                        given element
- *     element            (I) element to find neighbours to
- *     pElem              (O) array with elements loaded starting at pElem[0]
- *     exclusive          (I) T - exclude the given element in vElem
- *                            F - include the given element in vElem
- *     return             number of elements added to pElem
- *
- *   template <typename Elem>
- *   int all_element_face_neighbours(Elem *const element,
- *                                   std::vector<MElement*> &vElem,
- *                                   const bool exclusive = true)
- *                     -- find all the elements connected to the faces of a
- *                        given element
- *     element            (I) element to find neighbours to
- *     vElem              (O) vector with elements loaded via push_back
- *     exclusive          (I) T - exclude the given element in vElem
- *                            F - include the given element in vElem
- *     return             number of elements added to vElem
- *
- *   template <typename Elem>
- *   int all_element_face_neighbours(Elem *const element, MElement **pElem
- *                                   const bool exclusive = true)
- *                     -- find all the elements connected to the faces of a
- *                        given element
- *     element            (I) element to find neighbours to
- *     pElem              (O) array with elements loaded starting at pElem[0]
- *     exclusive          (I) T - exclude the given element in vElem
- *                            F - include the given element in vElem
- *     return             number of elements added to pElem
- *
- *   int element_vertex_neighbours(const MElement *const element,
- *                                 MVertex *const vertex,
- *                                 std::vector<MElement*> &vElem)
- *                     -- find all the elements connected to a vertex of a given
- *                        element.  This is the same as routine
- *                        vertex_neighbours except that the given element is
- *                        excluded from the result.
- *     element            (I) element to exclude
- *     vertex             (I) vertex to find element neighbours to
- *     vElem              (O) vector with elements loaded via push_back
- *     return             number of elements added to vElem
- *
- *   int element_vertex_neighbours(const MElement *const element,
- *                                 MVertex *const vertex, MElement **pElem)
- *                     -- find all the elements connected to a vertex of a given
- *                        element.  This is the same as routine
- *                        vertex_neighbours except that the given element is
- *                        excluded from the result.
- *     element            (I) element to exclude
- *     vertex             (I) vertex to find element neighbours to
- *     pElem              (O) array with elements loaded starting at pElem[0]
- *     return             number of elements added to pElem
- *
- *   int element_edge_neighbours(const MElement *const element,
- *                               const MEdge &edge,
- *                               std::vector<MElement*> &vElem)
- *                     -- find all the elements connected to an edge of a given
- *                        element.  This is the same as routine
- *                        edge_neighbours except that the given element is
- *                        excluded from the result.
- *     element            (I) element to exclude
- *     edge               (I) edge to find element neighbours to
- *     vElem              (O) vector with elements loaded via push_back
- *     return             number of elements added to vElem
- *
- *   int element_edge_neighbours(const MElement *const element,
- *                               const MEdge &edge, MElement **pElem)
- *                     -- find all the elements connected to an edge of a given
- *                        element.  This is the same as routine
- *                        edge_neighbours except that the given element is
- *                        excluded from the result.
- *     element            (I) element to exclude
- *     edge               (I) edge to find element neighbours to
- *     pElem              (O) array with elements loaded starting at pElem[0]
- *     return             number of elements added to pElem
- *
- *   MElement *element_face_neighbour(const MElement *const element,
- *                                    const MFace &face)
- *                     -- find the element connected to the face of a given
- *                        element.  Routine element_neighbour could also be
- *                        used.
- *     element            (I) element to find neighbour to
- *     face               (I) face to find element neighbour across
- *     return             the neighbour element or 0 if none found
- *
- *   MElement *element_neighbour(MElement *const element,
- *                               const int iPolytope)
- *                     -- find the neighbour element across the d-1 bounding
- *                        polytope.  E.g., for a triangle, find an element
- *                        across an edge.
- *     element            (I) element to find neighbour to
- *     iPolytope          (I) which polytope to find.  E.g., edge 1
- *     return             the neighbour element or 0 if none found
- *
- * Notes
- * =====
- *
- *   - Many of the query routines that have an element argument are specialized.
- *     If a specific element pointer is given, static (compile-time)
- *     polymorphism is used.  If a base element pointer is given, dynamic
- *     (run-time) polymorphism is used.  This is mostly experimental and may not
- *     have any tangible benefit.
- *
- * Example
- * =======
- *
- *------------------------------------------------------------------------------
-
-  // Examples for adding elements from containers (the type of container doesn't
-  // matter.
-
-  MNeighbour meshNeighbour;
-
-  std::list<GEntity*> faceEnt1;  // Assume known to be GFace type
-  meshNeighbour.add_elements_in_entities<GFace>(faceEnt1.begin(),
-                                                faceEnt1.end());
-
-  std::vector<GFace*> faceEnt2;
-  meshNeighbour.add_elements_in_entities(faceEnt2.begin(), faceEnt2.end());
-  meshNeighbour.add_elements_in_entity(faceEnt2.begin());
-
-  GEntity *faceEnt3[2];  // Assume known to be GFace type
-  meshNeighbour.add_elements_in_entities<GFace>(faceEnt3, 2);
-
-  GFace *faceEnt4[10];
-  meshNeighbour.add_elements_in_entities(faceEnt4, 10);
-
-  std::list<MElement*> elem1;
-  // Run time introspection.  Elements can be mixed.
-  meshNeighbour.add_elements(elem1.begin(), elem1.end());
-
-  std::vector<MTriangle*> elem2;
-  // Compile time introspection.
-  meshNeighbour.add_elements(elem2.begin(), elem2.end());
-
-  MElement *elem3[5];
-  // Run time introspection.  Elements can be mixed.
-  meshNeighbour.add_elements(elem3, 5);
-
-  MQuadrangle *elem4[6];
-  // Compile time introspection.
-  meshNeighbour.add_elements(elem4, 6);
-
- *------------------------------------------------------------------------------
- *
- ******************************************************************************/
-
-class MNeighbour
-{
-
-
-/*==============================================================================
- * Class scope types
- *============================================================================*/
-
- private:
-  typedef std::set<MElement*, std::less<MElement*> > ElemSet;
-  typedef ElemSet::const_iterator ElemSetConstIterator;
-
-
-/*==============================================================================
- * Iterators
- *============================================================================*/
-
- public:
-
-  // Iterators over polytopes
-  typedef PolytopeIterator<MVertex*> Vertex_const_iterator;
-  typedef PolytopeIterator<MEdge> Edge_const_iterator;
-  typedef PolytopeIterator<MFace> Face_const_iterator;
-
-  // Begin and end points for these iterators
-  Vertex_const_iterator vertex_begin()
-  {
-    return Vertex_const_iterator(vertNRange.begin()); 
-  }
-  Vertex_const_iterator vertex_end()
-  {
-    return Vertex_const_iterator(vertNRange.end()); 
-  }
-  Edge_const_iterator edge_begin()
-  {
-    return Edge_const_iterator(edgeNRange.begin()); 
-  }
-  Edge_const_iterator edge_end()
-  {
-    return Edge_const_iterator(edgeNRange.end()); 
-  }
-  Face_const_iterator face_begin()
-  {
-    return Face_const_iterator(faceNRange.begin()); 
-  }
-  Face_const_iterator face_end()
-  {
-    return Face_const_iterator(faceNRange.end()); 
-  }
-
-
-/*==============================================================================
- * Template meta-programming classes
- *============================================================================*/
-
- private:
-
-  template <typename Elem, typename Polytope, int NP>
-  struct FindNeighbours;
-  template <typename Ent, int N = EntTr<Ent>::numElemTypes>
-  struct FindNeighboursInEntity;
-
-
-/*==============================================================================
- * Public member functions and data
- *============================================================================*/
-
- public:
-
-//--Default constructor.
-
-  // The constructor cannot really do anything because the initialization
-  // functions (find_neighbours_*) cannot perform template argument deduction.
-  MNeighbour() : maxVertNeighbours(0), maxEdgeNeighbours(0),
-    maxFaceNeighbours(0) { }
-
-/*--------------------------------------------------------------------*
- * Elements added from entities.
- *--------------------------------------------------------------------*/
-
-//--Compute database of neighbour elements.  The entity pointers may be of a
-//--base or specific type and the call must state the specific type.
-
-  template <typename Ent, typename EntIter>
-  void add_elements_in_entities(EntIter begin, EntIter end) {
-
-    // Check the Ent and EntIter types
-    // NOTE:  If the compiler sent you here, you're using an invalid entity as a
-    // template parameter and/or and invalid iterator.  See struct
-    // ValidEntityIterator above for valid types
-    typedef typename ValidSpecEntityIterator<Ent,
-      typename EntIter::value_type>::Type Check;
-
-    // Find the neighbours of each vertex, edge, and face
-    for(EntIter itEnt = begin; itEnt != end; ++itEnt) {
-      Ent *entity = static_cast<Ent*>(*itEnt);
-      FindNeighboursInEntity<Ent>::eval(entity, vertNRange, edgeNRange,
-                                        faceNRange, neighbours);
-    }
-    maxVertNeighbours = -1;
-    maxEdgeNeighbours = -1;
-    maxFaceNeighbours = -1;
-  }
-
-//--Same as the routine above except the arguments are pointers instead of
-//--iterators.  Hence, value_type does not exist.
-
-  template <typename Ent, typename EntP>
-  void add_elements_in_entities(EntP **pEnt, unsigned int n) {
-
-    // Check the Ent and EntP types
-    // NOTE:  If the compiler sent you here, you're using an invalid entity as a
-    // template parameter and/or invalid pointers. See struct
-    // ValidEntityIterator above for valid types
-    typedef typename ValidSpecEntityIterator<Ent, EntP*>::Type Check;
-
-    // Find the neighbours of each vertex, edge, and face
-    while(n--) {
-      Ent *entity = static_cast<Ent*>(*pEnt++);
-      FindNeighboursInEntity<Ent>::eval(entity, vertNRange, edgeNRange,
-                                        faceNRange, neighbours);
-    }
-    maxVertNeighbours = -1;
-    maxEdgeNeighbours = -1;
-    maxFaceNeighbours = -1;
-  }
-
-//--Compute database of neighbour elements.  The entity pointers must be of a
-//--specific type (not a base class).  The type of entity is determined from
-//--Iterator::value_type and does not need to be stated.
-
-  template <typename EntIter>
-  void add_elements_in_entities(EntIter begin, EntIter end) {
-
-    // Check the value type in EntIter (it must not be GEntity*)
-    // NOTE:  If the compiler sent you here, you're using an invalid entity as a
-    // template parameter and/or and invalid iterator for the constructor.  See
-    // struct ValidEntityIterator above for valid types
-    typedef typename ValidEntityIterator<typename EntIter::value_type>::
-      Type Ent;
-
-    // Find the neighbours of each vertex, edge, and face
-    for(EntIter itEnt = begin; itEnt != end; ++itEnt) {
-      FindNeighboursInEntity<Ent>::eval(*itEnt, vertNRange, edgeNRange,
-                                        faceNRange, neighbours);
-    }
-    maxVertNeighbours = -1;
-    maxEdgeNeighbours = -1;
-    maxFaceNeighbours = -1;
-  }
-
-//--Same as the routine above except the arguments are pointers instead of
-//--iterators.  Hence, value_type does not exist.
-
-  template <typename Ent>
-  void add_elements_in_entities(Ent **pEnt, unsigned int n) {
-
-    // Check the type Ent (it must not be GEntity)
-    // NOTE:  If the compiler sent you here, you're using an invalid entity as a
-    // pointer See struct ValidEntityIterator above for valid types
-    typedef typename ValidEntityIterator<Ent*>::Type Check;
-
-    // Find the neighbours of each vertex, edge, and face
-    while(n--) {
-      FindNeighboursInEntity<Ent>::eval(*pEnt++, vertNRange, edgeNRange,
-                                        faceNRange, neighbours);
-    }
-    maxVertNeighbours = -1;
-    maxEdgeNeighbours = -1;
-    maxFaceNeighbours = -1;
-  }
-
-//--Find neighbours from the elements attached to a pointer to a single entity.
-//--The type of entity does not need to be stated.
-
-  template <typename EntPtr>
-  void add_elements_in_entity(EntPtr entity) {
-
-    // Check the type of EntPtr (it must not be GEntity*)
-    // NOTE:  If the compiler sent you here, you're using an invalid entity as a
-    // pointer See struct ValidEntityIterator above for valid types
-    typedef typename ValidEntityIterator<EntPtr>::Type Ent;
-
-    // Find the neighbours of each vertex, edge, and face
-    FindNeighboursInEntity<Ent>::eval(entity, vertNRange, edgeNRange,
-                                      faceNRange, neighbours);
-    maxVertNeighbours = -1;
-    maxEdgeNeighbours = -1;
-    maxFaceNeighbours = -1;
-  }
-
-/*--------------------------------------------------------------------*
- * Elements added directly
- *--------------------------------------------------------------------*/
-
-//--Add elements from a container using iterators
-
-  template <typename ElemIter>
-  void add_elements(ElemIter begin, ElemIter end) {
-
-    // Strip the pointer to the element type and convert it to first order
-    typedef typename ElemBaseOrder<
-      typename StripPointer<
-        typename ElemIter::value_type
-      >::Type
-    >::Order1 Elem_1o;
-
-    // Check for a valid type of element
-    // NOTE:  If the compiler sent you here, you're using an invalid iterator.
-    // The dereferenced object needs to be a pointer to an element.
-    typedef typename ValidElement<Elem_1o>::Type Check;
-
-    add_elements1<Elem_1o, ElemIter>(begin, end);
-  }
-
-//--Add elements from an array of pointers (value_type does not exist)
-
-  template <typename Elem>
-  void add_elements(Elem **pElem, const unsigned int n) {
-
-    // Convert the element type to first order
-    typedef typename ElemBaseOrder<Elem>::Order1 Elem_1o;
-
-    // Check for a valid type of element
-    // NOTE:  If the compiler sent you here, you're using an invalid pointer.
-    // The dereferenced object needs to be a pointer to an element.
-    typedef typename ValidElement<Elem_1o>::Type Check;
-
-    add_elements1<Elem_1o, Elem**>(pElem, pElem+n);
-  }
-
-/*--------------------------------------------------------------------*
- * Reset the database
- *--------------------------------------------------------------------*/
-
-  void clear() 
-  {
-    vertNRange.clear();
-    edgeNRange.clear();
-    faceNRange.clear();
-    neighbours.clear();
-    maxVertNeighbours = 0;
-    maxEdgeNeighbours = 0;
-    maxFaceNeighbours = 0;
-  }
-
-/*--------------------------------------------------------------------*
- * Query the database
- *--------------------------------------------------------------------*/
-
-//--Get the max number of elements sharing a vertex, edge, or face
-
-  int max_vertex_neighbours()
-  {
-     if ( maxVertNeighbours < 0 ) compute_max_vert_neighbours();
-     return maxVertNeighbours;
-  }
-  int max_edge_neighbours()
-  {
-     if ( maxEdgeNeighbours < 0 ) compute_max_edge_neighbours();
-     return maxEdgeNeighbours;
-  }
-  int max_face_neighbours()
-  {
-     if ( maxFaceNeighbours < 0 ) compute_max_face_neighbours();
-     return maxFaceNeighbours;
-  }
-
-//--Return elements that share a vertex
-
-  // From a vertex iterator
-  int vertex_neighbours(Vertex_const_iterator &itVert,
-                        std::vector<MElement*> &vElem) const
-  {
-    NeighboursConstIterator itElem = itVert.begin_neighbours();
-    for(int n = itVert.num_neighbours(); n--;) vElem.push_back(*itElem++);
-    return itVert.num_neighbours();
-  }
-
-  int vertex_neighbours(Vertex_const_iterator &itVert, MElement **pElem) const
-  {
-    NeighboursConstIterator itElem = itVert.begin_neighbours();
-    for(int n = itVert.num_neighbours(); n--;) *pElem++ = *itElem++;
-    return itVert.num_neighbours();
-  }
-
-  // From a vertex
-  int vertex_neighbours(MVertex *const vertex,
-                        std::vector<MElement*> &vElem) const
-  {
-    VertNRangeConstIterator itVert = vertNRange.find(vertex);
-    if(itVert == vertNRange.end()) return 0;  // Vertex not found
-    NeighboursConstIterator itElem = itVert->second.begin;
-    for(int n = itVert->second.num; n--;) vElem.push_back(*itElem++);
-    return itVert->second.num;
-  }
-
-  int vertex_neighbours(MVertex *const vertex, MElement **pElem) const
-  {
-    VertNRangeConstIterator itVert = vertNRange.find(vertex);
-    if(itVert == vertNRange.end()) return 0;  // Vertex not found
-    NeighboursConstIterator itElem = itVert->second.begin;
-    for(int n = itVert->second.num; n--;) *pElem++ = *itElem++;
-    return itVert->second.num;
-  }
-
-//--Return elements that share an edge
-
-  // From an edge iterator
-  int edge_neighbours(Edge_const_iterator &itEdge,
-                      std::vector<MElement*> &vElem) const 
-  {
-    NeighboursConstIterator itElem = itEdge.begin_neighbours();
-    for(int n = itEdge.num_neighbours(); n--;)
-       vElem.push_back(*itElem++);
-    return itEdge.num_neighbours();
-  }
-
-  int edge_neighbours(Edge_const_iterator &itEdge, MElement **pElem) const
-  {
-    NeighboursConstIterator itElem = itEdge.begin_neighbours();
-    for(int n = itEdge.num_neighbours(); n--;)
-       *pElem++ = *itElem++;
-    return itEdge.num_neighbours();
-  }
-
-  // From an edge
-  int edge_neighbours(const MEdge &edge, std::vector<MElement*> &vElem) const
-  {
-    EdgeNRangeConstIterator itEdge = edgeNRange.find(edge);
-    if(itEdge == edgeNRange.end()) return 0;  // Edge not found
-    NeighboursConstIterator itElem = itEdge->second.begin;
-    for(int n = itEdge->second.num; n--;) vElem.push_back(*itElem++);
-    return itEdge->second.num;
-  }
-
-  int edge_neighbours(const MEdge &edge, MElement **pElem) const
-  {
-    EdgeNRangeConstIterator itEdge = edgeNRange.find(edge);
-    if(itEdge == edgeNRange.end()) return 0;  // Edge not found
-    NeighboursConstIterator itElem = itEdge->second.begin;
-    for(int n = itEdge->second.num; n--;) *pElem++ = *itElem++;
-    return itEdge->second.num;
-  }
-
-//--Return elements that share a face
-
-  // From a face iterator
-  int face_neighbours(Face_const_iterator &itFace,
-                      std::vector<MElement*> &vElem) const 
-  {
-    NeighboursConstIterator itElem = itFace.begin_neighbours();
-    for(int n = itFace.num_neighbours(); n--;)
-       vElem.push_back(*itElem++);
-    return itFace.num_neighbours();
-  }
-
-  int face_neighbours(Face_const_iterator &itFace, MElement **pElem) const
-  {
-    NeighboursConstIterator itElem = itFace.begin_neighbours();
-    for(int n = itFace.num_neighbours(); n--;)
-       *pElem++ = *itElem++;
-    return itFace.num_neighbours();
-  }
-
-  // From a face
-  int face_neighbours(const MFace &face, std::vector<MElement*> &vElem) const
-  {
-    FaceNRangeConstIterator itFace = faceNRange.find(face);
-    if(itFace == faceNRange.end()) return 0;  // Face not found
-    NeighboursConstIterator itElem = itFace->second.begin;
-    for(int n = itFace->second.num; n--;) vElem.push_back(*itElem++);
-    return itFace->second.num;
-  }
-
-  int face_neighbours(const MFace &face, MElement **pElem) const
-  {
-    FaceNRangeConstIterator itFace = faceNRange.find(face);
-    if(itFace == faceNRange.end()) return 0;  // Face not found
-    NeighboursConstIterator itElem = itFace->second.begin;
-    for(int n = itFace->second.num; n--;) *pElem++ = *itElem++;
-    return itFace->second.num;
-  }
-
-//--Return the number of elements that share a vertex
-
-  int vertex_num_neighbours(MVertex *const vertex) const
-  {
-    VertNRangeConstIterator itVert = vertNRange.find(vertex);
-    if(itVert == vertNRange.end()) return 0;  // Vertex not found
-    return itVert->second.num;
-  }
-
-//--Return the number of elements that share an edge
-
-  int edge_num_neighbours(const MEdge& edge) const
-  {
-    EdgeNRangeConstIterator itEdge = edgeNRange.find(edge);
-    if(itEdge == edgeNRange.end()) return 0;  // Edge not found
-    return itEdge->second.num;
-  }
-
-//--Return the number of elements that share a face
-
-  int face_num_neighbours(const MFace& face) const
-  {
-    FaceNRangeConstIterator itFace = faceNRange.find(face);
-    if(itFace == faceNRange.end()) return 0;  // Face not found
-    return itFace->second.num;
-  }
-
-//--Return all the elements connected to a given element (Note: these routines
-//--have specializations - they are defined following the class)
-
-  template <typename Elem>
-  int all_element_neighbours(Elem *const element, std::vector<MElement*> &vElem,
-                             const bool exclusive = true)
-  {
-    typedef typename ElemBaseOrder<Elem>::Order1 Elem_1o;
-    element_neighbours1<Elem, MVertex*>(static_cast<Elem_1o*>(element),
-                                        vertNRange, ElemTr<Elem_1o>::numVertex);
-    if(exclusive) elemSet.erase(element);
-    const ElemSetConstIterator itElemEnd = elemSet.end();
-    for(ElemSetConstIterator itElem = elemSet.begin(); itElem != itElemEnd;
-        ++itElem) vElem.push_back(*itElem);
-    const int nElem = elemSet.size();
-    elemSet.clear();
-    return nElem;
-  }
-
-  template <typename Elem>
-  int all_element_neighbours(Elem *const element, MElement **pElem,
-                             const bool exclusive = true)
-  {
-    typedef typename ElemBaseOrder<Elem>::Order1 Elem_1o;
-    element_neighbours1<Elem, MVertex*>(static_cast<Elem_1o*>(element),
-                                        vertNRange, ElemTr<Elem_1o>::numVertex);
-    if(exclusive) elemSet.erase(element);
-    const ElemSetConstIterator itElemEnd = elemSet.end();
-    for(ElemSetConstIterator itElem = elemSet.begin(); itElem != itElemEnd;
-        ++itElem) *pElem++ = *itElem;
-    const int nElem = elemSet.size();
-    elemSet.clear();
-    return nElem;
-  }
-
-//--Return all the elements connected to the edges of a given element (Note:
-//--these routines have specializations - they are defined following the class)
-
-  template <typename Elem>
-  int all_element_edge_neighbours(Elem *const element,
-                                  std::vector<MElement*> &vElem,
-                                  const bool exclusive = true)
-  {
-    typedef typename ElemBaseOrder<Elem>::Order1 Elem_1o;
-    element_neighbours1<Elem, MEdge>(static_cast<Elem_1o*>(element),
-                                     edgeNRange, ElemTr<Elem_1o>::numEdge);
-    if(exclusive) elemSet.erase(element);
-    const ElemSetConstIterator itElemEnd = elemSet.end();
-    for(ElemSetConstIterator itElem = elemSet.begin(); itElem != itElemEnd;
-        ++itElem) vElem.push_back(*itElem);
-    const int nElem = elemSet.size();
-    elemSet.clear();
-    return nElem;
-  }
-
-  template <typename Elem>
-  int all_element_edge_neighbours(Elem *const element, MElement **pElem,
-                                  const bool exclusive = true)
-  {
-    typedef typename ElemBaseOrder<Elem>::Order1 Elem_1o;
-    element_neighbours1<Elem, MEdge>(static_cast<Elem_1o*>(element),
-                                     edgeNRange, ElemTr<Elem_1o>::numEdge);
-    if(exclusive) elemSet.erase(element);
-    const ElemSetConstIterator itElemEnd = elemSet.end();
-    for(ElemSetConstIterator itElem = elemSet.begin(); itElem != itElemEnd;
-        ++itElem) *pElem++ = *itElem;
-    const int nElem = elemSet.size();
-    elemSet.clear();
-    return nElem;
-  }
-
-//--Return all the elements connected to the faces of a given element (Note:
-//--these routines have specializations - they are defined following the class)
-
-  template <typename Elem>
-  int all_element_face_neighbours(Elem *const element,
-                                  std::vector<MElement*> &vElem,
-                                  const bool exclusive = true)
-  {
-    typedef typename ElemBaseOrder<Elem>::Order1 Elem_1o;
-    element_neighbours1<Elem, MFace>(static_cast<Elem_1o*>(element),
-                                     faceNRange, ElemTr<Elem_1o>::numFace);
-    if(exclusive) elemSet.erase(element);
-    const ElemSetConstIterator itElemEnd = elemSet.end();
-    for(ElemSetConstIterator itElem = elemSet.begin(); itElem != itElemEnd;
-        ++itElem) vElem.push_back(*itElem);
-    const int nElem = elemSet.size();
-    elemSet.clear();
-    return nElem;
-  }
-
-  template <typename Elem>
-  int all_element_face_neighbours(Elem *const element, MElement **pElem,
-                                  const bool exclusive = true)
-  {
-    typedef typename ElemBaseOrder<Elem>::Order1 Elem_1o;
-    element_neighbours1<Elem, MFace>(static_cast<Elem_1o*>(element),
-                                     faceNRange, ElemTr<Elem_1o>::numFace);
-    if(exclusive) elemSet.erase(element);
-    const ElemSetConstIterator itElemEnd = elemSet.end();
-    for(ElemSetConstIterator itElem = elemSet.begin(); itElem != itElemEnd;
-        ++itElem) *pElem++ = *itElem;
-    const int nElem = elemSet.size();
-    elemSet.clear();
-    return nElem;
-  }
-
-//--Return elements connected to a specific vertex of a given element (Note:
-//--these routines are the same as vertex_neighbours() except that the given
-//--element is excluded)
-
-  int element_vertex_neighbours(const MElement *const element,
-                                MVertex *const vertex,
-                                std::vector<MElement*> &vElem) const
-  {
-    VertNRangeConstIterator itVert = vertNRange.find(vertex);
-    if(itVert == vertNRange.end()) return 0;  // Vertex not found
-    NeighboursConstIterator itElem = itVert->second.begin;
-    int n = itVert->second.num;
-    while(*itElem != element) {
-      vElem.push_back(*itElem++);
-      --n;
-    }
-    ++itElem;
-    for(--n; n--;) vElem.push_back(*itElem++);
-    return itVert->second.num - 1;
-  }
-
-  int element_vertex_neighbours(const MElement *const element,
-                                MVertex *const vertex, MElement **pElem) const
-  {
-    VertNRangeConstIterator itVert = vertNRange.find(vertex);
-    if(itVert == vertNRange.end()) return 0;  // Vertex not found
-    NeighboursConstIterator itElem = itVert->second.begin;
-    int n = itVert->second.num;
-    while(*itElem != element) {
-      *pElem++ = *itElem++;
-      --n;
-    }
-    ++itElem;
-    for(--n; n--;) *pElem++ = *itElem++;
-    return itVert->second.num - 1;
-  }
-
-//--Return elements connected to a specific edge of a given element (Note: these
-//--routines are the same as edge_neighbours except that the given element is
-//--excluded)
-
-  int element_edge_neighbours(const MElement *const element, const MEdge &edge,
-                              std::vector<MElement*> &vElem) const
-  {
-    EdgeNRangeConstIterator itEdge = edgeNRange.find(edge);
-    if(itEdge == edgeNRange.end()) return 0;  // Edge not found
-    NeighboursConstIterator itElem = itEdge->second.begin;
-    int n = itEdge->second.num;
-    while(*itElem != element) {
-      vElem.push_back(*itElem++);
-      --n;
-    }
-    ++itElem;
-    for(--n; n--;) vElem.push_back(*itElem++);
-    return itEdge->second.num - 1;
-  }
-
-  int element_edge_neighbours(const MElement *const element, const MEdge &edge,
-                              MElement **pElem) const
-  {
-    EdgeNRangeConstIterator itEdge = edgeNRange.find(edge);
-    if(itEdge == edgeNRange.end()) return 0;  // Edge not found
-    NeighboursConstIterator itElem = itEdge->second.begin;
-    int n = itEdge->second.num;
-    while(*itElem != element) {
-      *pElem++ = *itElem++;
-      --n;
-    }
-    ++itElem;
-    for(--n; n--;) *pElem++ = *itElem++;
-    return itEdge->second.num - 1;
-  }
-
-//--Return the element connected to a specific face of a given element (Note:
-//--this routine is similar to element_neighbour except less general)
-
-  MElement *element_face_neighbour(const MElement *const element,
-                                   const MFace &face) const
-  {
-    FaceNRangeConstIterator itFace = faceNRange.find(face);
-    if(itFace == faceNRange.end()) return 0;  // Face not found
-    NeighboursConstIterator itElem = itFace->second.begin;
-    if(*itElem != element) return *itElem;
-    if(itFace->second.num > 1) {
-      ++itElem;
-      if(*itElem != element) return *itElem;
-    }
-    return 0;
-  }
-
-//--Return the element neighbour across a d-1 polytope
-
-  MElement *element_neighbour(MElement *const element,
-                              const int iPolytope) const
-  {
-    int num;
-    NeighboursConstIterator itElem;
-    switch(element->getTypeForMSH()) {
-    case MSH_LIN_2:
-    case MSH_LIN_3:
-      {
-        MVertex *vertex = element->getVertex(iPolytope);
-        VertNRangeConstIterator itVert = vertNRange.find(vertex);
-        if(itVert == vertNRange.end()) return 0;
-        num = itVert->second.num;
-        itElem = itVert->second.begin;
-      }
-      break;
-    case MSH_TRI_3:
-    case MSH_TRI_6:
-    case MSH_QUA_4:
-    case MSH_QUA_8:
-    case MSH_QUA_9:
-      {
-        MEdge edge = element->getEdge(iPolytope);
-        EdgeNRangeConstIterator itEdge = edgeNRange.find(edge);
-        if(itEdge == edgeNRange.end()) return 0;
-        num = itEdge->second.num;
-        itElem = itEdge->second.begin;
-      }
-      break;
-    case MSH_TET_4:
-    case MSH_TET_10:
-    case MSH_HEX_8:
-    case MSH_HEX_20:
-    case MSH_HEX_27:
-    case MSH_PRI_6:
-    case MSH_PRI_15:
-    case MSH_PRI_18:
-    case MSH_PYR_5:
-    case MSH_PYR_13:
-    case MSH_PYR_14: {
-        MFace face = element->getFace(iPolytope);
-        FaceNRangeConstIterator itFace = faceNRange.find(face);
-        if(itFace == faceNRange.end()) return 0;
-        num = itFace->second.num;
-        itElem = itFace->second.begin;
-      }
-      break;
-    }
-    if(*itElem != element) return *itElem;
-    if(num > 1) {
-      ++itElem;
-      if(*itElem != element) return  *itElem;
-    }
-    return 0;
-  }
-
-
-/*==============================================================================
- * Private member functions and data
- *============================================================================*/
-
-private:
-
-//--Data members
-
-  VertNRange vertNRange;                // Map of Vert. to indexes in neighbours
-  EdgeNRange edgeNRange;                // Map of edge to indexes in neighbours
-  FaceNRange faceNRange;                // Map of face to indexes in neighbours
-  Neighbours neighbours;                // List of element neighbours
-  ElemSet elemSet;                      // Buffer for finding unique elements
-  int maxVertNeighbours;
-  int maxEdgeNeighbours;
-  int maxFaceNeighbours;
-
-//--Working routine to determine the neighbours from elements (Note: Ideally
-//--this class would use a specialization for MElement but partial
-//--specialization is not currently permitted for member functions)
-
-  template <typename Elem, typename ElemIter>
-  void add_elements1(ElemIter begin, ElemIter end)
-  {
-    // Test if the type of element is known
-    if(ElemTr<Elem>::numVertex) {  // 0 only for MElement
-      // Find the neighbours of each vertex, edge, and face
-      for(ElemIter itElem = begin; itElem != end; ++itElem)
-        find_polytope_neighbours<Elem>(*itElem);
-    }
-    else {
-      // We only have an iterator to an MElement* and have to check the type in
-      // run-time.
-      for(ElemIter itElem = begin; itElem != end; ++itElem) {
-        switch((*itElem)->getTypeForMSH()) {
-        case MSH_LIN_2:
-        case MSH_LIN_3:
-          find_polytope_neighbours<MLine>(*itElem);
-          break;
-        case MSH_TRI_3:
-        case MSH_TRI_6:
-          find_polytope_neighbours<MTriangle>(*itElem);
-          break;
-        case MSH_QUA_4:
-        case MSH_QUA_8:
-        case MSH_QUA_9:
-          find_polytope_neighbours<MQuadrangle>(*itElem);
-          break;
-        case MSH_TET_4:
-        case MSH_TET_10:
-          find_polytope_neighbours<MTetrahedron>(*itElem);
-          break;
-        case MSH_HEX_8:
-        case MSH_HEX_20:
-        case MSH_HEX_27:
-          find_polytope_neighbours<MHexahedron>(*itElem);
-          break;
-        case MSH_PRI_6:
-        case MSH_PRI_15:
-        case MSH_PRI_18:
-          find_polytope_neighbours<MPrism>(*itElem);
-          break;
-        case MSH_PYR_5:
-        case MSH_PYR_13:
-        case MSH_PYR_14:
-          find_polytope_neighbours<MPyramid>(*itElem);
-          break;
-        }
-      }
-    }
-    maxVertNeighbours = -1;
-    maxEdgeNeighbours = -1;
-    maxFaceNeighbours = -1;
-  }
-
-//--Find the neighbours sharing each polytope
-
-  template <typename Elem>
-  void find_polytope_neighbours(MElement *element)
-  {
-    // Find neighbours for vertices
-    FindNeighbours<Elem, MVertex*, ElemTr<Elem>::numVertex>::
-      eval(static_cast<Elem*>(element), vertNRange, neighbours);
-    // Find neighbours for edges
-    FindNeighbours<Elem, MEdge, ElemTr<Elem>::numEdge>::
-      eval(static_cast<Elem*>(element), edgeNRange, neighbours);
-    // Find neighbours for faces
-    FindNeighbours<Elem, MFace, ElemTr<Elem>::numFace>::
-      eval(static_cast<Elem*>(element), faceNRange, neighbours);
-  }
-
-//--Compute max neighbours
-
-  void compute_max_vert_neighbours()
-  {
-    maxVertNeighbours = 0;
-    for(VertNRangeConstIterator itVert = vertNRange.begin();
-        itVert != vertNRange.end(); ++itVert)
-      maxVertNeighbours = std::max(maxVertNeighbours, itVert->second.num);
-  }
-  void compute_max_edge_neighbours()
-  {
-    maxEdgeNeighbours = 0;
-    for(EdgeNRangeConstIterator itEdge = edgeNRange.begin();
-        itEdge != edgeNRange.end(); ++itEdge)
-      maxEdgeNeighbours = std::max(maxEdgeNeighbours, itEdge->second.num);
-  }
-  void compute_max_face_neighbours()
-  {
-    maxFaceNeighbours = 0;
-    for(FaceNRangeConstIterator itFace = faceNRange.begin();
-        itFace != faceNRange.end(); ++itFace)
-      maxFaceNeighbours = std::max(maxFaceNeighbours, itFace->second.num);
-  }
-
-//--Work routine for finding all the element neighbours
-
-  template <typename Elem, typename Polytope>
-  void element_neighbours1(Elem *const element,
-                           typename PolytopeTr<Polytope>::PolytopeNRange
-                           &polytopeNRange, const int nVPE)
-  {
-    for(int iVPE = 0; iVPE != nVPE; ++iVPE) {
-      Polytope polytope = PolytopeTr<Polytope>::getPolytope(element, iVPE);
-      typename PolytopeTr<Polytope>::PNRConstIterator itPoly =
-        polytopeNRange.find(polytope);
-      if(itPoly != polytopeNRange.end()) {
-        NeighboursConstIterator itElem = itPoly->second.begin;
-        for(int n = itPoly->second.num; n--;)
-          elemSet.insert(*itElem++);
-      }
-    }
-  }
-
-};
-
-
-/*==============================================================================
- * Template meta-programming classes
- *============================================================================*/
-
-/*--------------------------------------------------------------------
- * This is the main working routine for finding the elements sharing
- * a lower-dimension polytope (i.e, the structures bounding the
- * element: vertex, edge, or face).  It is templated for any type of
- * polytope and uses template meta-programming to unroll the loop over
- * the number of that type of polytope in the element.
- *--------------------------------------------------------------------*/
-
-//--Entry point
-
-template <typename Elem,                // Elem is the type of element
-          typename Polytope,            // Polytope is a lower dimensional
-                                        // bounding object (MVertex*, MEdge, or
-                                        // MFace) of the element.  We want to
-                                        // find the elements sharing a given
-                                        // Polytope.
-          int NP>                       // NP is an index through the number
-                                        // of Polytope in the element (e.g.,
-                                        // number of vertices)
-struct MNeighbour::FindNeighbours
-{
-  typedef typename PolytopeTr<Polytope>::PolytopeNRange PolytopeNRange;
-  static void eval(Elem *element, PolytopeNRange &polytopeNRange,
-                   Neighbours &neighbours)
-  {
-    Polytope polytope = PolytopeTr<Polytope>::template
-      getPolytope<Elem>(element, NP-1);
-    Range_t &range = polytopeNRange[polytope];
-    if(range.num == 0) range.begin = neighbours.end();
-    range.begin = neighbours.insert(range.begin, element);
-    ++range.num;
-    FindNeighbours<Elem, Polytope, NP-1>::eval(element, polytopeNRange,
-                                               neighbours);
-  }
-};
-
-//--Terminate loop when NP == 1
-
-template <typename Elem, typename Polytope>
-struct MNeighbour::FindNeighbours<Elem, Polytope, 1>
-{
-  typedef typename PolytopeTr<Polytope>::PolytopeNRange PolytopeNRange;
-  static void eval(Elem *element, PolytopeNRange &polytopeNRange,
-                   Neighbours &neighbours)
-  {
-    Polytope polytope = PolytopeTr<Polytope>::template
-      getPolytope<Elem>(element, 0);
-    Range_t &range = polytopeNRange[polytope];
-    if(range.num == 0) range.begin = neighbours.end();
-    range.begin = neighbours.insert(range.begin, element);
-    ++range.num;
-  }
-};
-
-//--Nothing to do if the dimension of this polytope is equal to or greater than
-//--that of the element (e.g., no faces (2D) in a 2D element).  This is
-//--indicated by NP = 0.
-
-template <typename Elem, typename Polytope>
-struct MNeighbour::FindNeighbours<Elem, Polytope, 0>
-{
-  typedef typename PolytopeTr<Polytope>::PolytopeNRange PolytopeNRange;
-  static void eval(Elem *element, PolytopeNRange &polytopeNRange,
-                   Neighbours &neighbours) { }
-};
-
-/*--------------------------------------------------------------------*
- * This class uses traits to determine the element types in an entity
- * and iterate over them.  A template meta-programming loop is used to
- * iterate over the types of elements.
- *--------------------------------------------------------------------*/
-
-//--Entry point
-
-template <typename Ent, int N>
-struct MNeighbour::FindNeighboursInEntity {
-  typedef typename EntElemTr<Ent, N>::Elem Elem;
-  static void eval(const Ent *const entity, VertNRange &vertNRange,
-                   EdgeNRange &edgeNRange, FaceNRange &faceNRange,
-                   Neighbours &neighbours)
-  {
-    for(typename EntElemTr<Ent, N>::ConstElementIterator itElem =
-          EntElemTr<Ent, N>::begin(entity);
-        itElem != EntElemTr<Ent, N>::end(entity); ++itElem) {
-      // Find neighbours for vertices
-      FindNeighbours<Elem, MVertex*, ElemTr<Elem>::numVertex>::
-        eval(*itElem, vertNRange, neighbours);
-      // Find neighbours for edges
-      FindNeighbours<Elem, MEdge, ElemTr<Elem>::numEdge>::
-        eval(*itElem, edgeNRange, neighbours);
-      // Find neighbours for faces
-      FindNeighbours<Elem, MFace, ElemTr<Elem>::numFace>::
-        eval(*itElem, faceNRange, neighbours);
-    }
-    FindNeighboursInEntity<Ent, N-1>::eval(entity, vertNRange, edgeNRange,
-                                           faceNRange, neighbours);
-  }
-};
-
-//--Terminate loop when N = 1
-
-template <typename Ent>
-struct MNeighbour::FindNeighboursInEntity<Ent, 1> {
-  typedef typename EntElemTr<Ent, 1>::Elem Elem;
-  static void eval(const Ent *const entity, VertNRange &vertNRange,
-                   EdgeNRange &edgeNRange, FaceNRange &faceNRange,
-                   Neighbours &neighbours)
-  {
-    for(typename EntElemTr<Ent, 1>::ConstElementIterator itElem =
-          EntElemTr<Ent, 1>::begin(entity);
-        itElem != EntElemTr<Ent, 1>::end(entity); ++itElem) {
-      // Find neighbours for vertices
-      FindNeighbours<Elem, MVertex*, ElemTr<Elem>::numVertex>::
-        eval(*itElem, vertNRange, neighbours);
-      // Find neighbours for edges
-      FindNeighbours<Elem, MEdge, ElemTr<Elem>::numEdge>::
-        eval(*itElem, edgeNRange, neighbours);
-      // Find neighbours for faces
-      FindNeighbours<Elem, MFace, ElemTr<Elem>::numFace>::
-        eval(*itElem, faceNRange, neighbours);
-    }
-  }
-};
-
-
-/*==============================================================================
- * Specializations for base elements - dynamic polymorphism
- *============================================================================*/
-
-//--Return all the elements connected to a given element
-
-template <>
-inline int MNeighbour::all_element_neighbours<MElement>
-(MElement *const element, std::vector<MElement*> &vElem, const bool exclusive)
-{
-  element_neighbours1<MElement, MVertex*>(element, vertNRange,
-                                          element->getNumPrimaryVertices());
-  if(exclusive) elemSet.erase(element);
-  const ElemSetConstIterator itElemEnd = elemSet.end();
-  for(ElemSetConstIterator itElem = elemSet.begin(); itElem != itElemEnd;
-      ++itElem) vElem.push_back(*itElem);
-  const int nElem = elemSet.size();
-  elemSet.clear();
-  return nElem;
-}
-
-template <>
-inline int MNeighbour::all_element_neighbours<MElement>
-(MElement *const element, MElement **pElem, const bool exclusive)
-{
-  element_neighbours1<MElement, MVertex*>(element, vertNRange,
-                                          element->getNumPrimaryVertices());
-  if(exclusive) elemSet.erase(element);
-  const ElemSetConstIterator itElemEnd = elemSet.end();
-  for(ElemSetConstIterator itElem = elemSet.begin(); itElem != itElemEnd;
-      ++itElem) *pElem++ = *itElem;
-  const int nElem = elemSet.size();
-  elemSet.clear();
-  return nElem;
-}
-
-//--Return all the elements connected to the edges of a given element
-
-template <>
-inline int MNeighbour::all_element_edge_neighbours<MElement>
-(MElement *const element, std::vector<MElement*> &vElem, const bool exclusive)
-{
-  element_neighbours1<MElement, MEdge>(element, edgeNRange,
-                                       element->getNumEdges());
-  if(exclusive) elemSet.erase(element);
-  const ElemSetConstIterator itElemEnd = elemSet.end();
-  for(ElemSetConstIterator itElem = elemSet.begin(); itElem != itElemEnd;
-      ++itElem) vElem.push_back(*itElem);
-  const int nElem = elemSet.size();
-  elemSet.clear();
-  return nElem;
-}
-
-template <>
-inline int MNeighbour::all_element_edge_neighbours<MElement>
-(MElement *const element, MElement **pElem, const bool exclusive)
-{
-  element_neighbours1<MElement, MEdge>(element, edgeNRange,
-                                       element->getNumEdges());
-  if(exclusive) elemSet.erase(element);
-  const ElemSetConstIterator itElemEnd = elemSet.end();
-  for(ElemSetConstIterator itElem = elemSet.begin(); itElem != itElemEnd;
-      ++itElem) *pElem++ = *itElem;
-  const int nElem = elemSet.size();
-  elemSet.clear();
-  return nElem;
-}
-
-//--Return all the elements connected to the faces of a given element
-
-template <>
-inline int MNeighbour::all_element_face_neighbours<MElement>
-(MElement *const element, std::vector<MElement*> &vElem, const bool exclusive)
-{
-  element_neighbours1<MElement, MFace>(element, faceNRange,
-                                       element->getNumFaces());
-  if(exclusive) elemSet.erase(element);
-  const ElemSetConstIterator itElemEnd = elemSet.end();
-  for(ElemSetConstIterator itElem = elemSet.begin(); itElem != itElemEnd;
-      ++itElem) vElem.push_back(*itElem);
-  const int nElem = elemSet.size();
-  elemSet.clear();
-  return nElem;
-}
-
-template <>
-inline int MNeighbour::all_element_face_neighbours<MElement>
-(MElement *const element, MElement **pElem, const bool exclusive)
-{
-  element_neighbours1<MElement, MFace>(element, faceNRange,
-                                       element->getNumFaces());
-  if(exclusive) elemSet.erase(element);
-  const ElemSetConstIterator itElemEnd = elemSet.end();
-  for(ElemSetConstIterator itElem = elemSet.begin(); itElem != itElemEnd;
-      ++itElem) *pElem++ = *itElem;
-  const int nElem = elemSet.size();
-  elemSet.clear();
-  return nElem;
-}
-
-#endif
diff --git a/Geo/MVertex.h b/Geo/MVertex.h
index 8c1a4a75c8..3d663fc027 100644
--- a/Geo/MVertex.h
+++ b/Geo/MVertex.h
@@ -74,7 +74,7 @@ class MVertex{
   inline double & x() { return _x; }
   inline double & y() { return _y; }
   inline double & z() { return _z; }
-  inline SPoint3 point() { return SPoint3(_x, _y, _z); }
+  inline SPoint3 point() const { return SPoint3(_x, _y, _z); }
 
   // get/set the parent entity
   inline GEntity* onWhat() const { return _ge; }
diff --git a/Geo/MZone.cpp b/Geo/MZone.cpp
new file mode 100644
index 0000000000..a0a9963e3a
--- /dev/null
+++ b/Geo/MZone.cpp
@@ -0,0 +1,396 @@
+// Gmsh - Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
+//
+// See the LICENSE.txt file for license information. Please report all
+// bugs and problems to <gmsh@geuz.org>.
+//
+// MZone.cpp - Copyright (C) 2008 S. Guzik, C. Geuzaine, J.-F. Remacle
+
+#include <iostream> // DBG
+#include "MZone.h"
+
+/*==============================================================================
+ * Forward declarations
+ *============================================================================*/
+
+template <typename Ent, int N = EntTr<Ent>::numElemTypes>
+struct ParseEntity;
+
+
+/*******************************************************************************
+ *
+ * Routines from class MZone
+ *
+ ******************************************************************************/
+
+
+/*==============================================================================
+ *
+ * Routine: add_elements_in_entities
+ *
+ * Purpose
+ * =======
+ *
+ *   Adds all (or only those of the right partition) elements in a container of
+ *   entities to the zone.
+ *
+ * I/O
+ * ===
+ *
+ *   <Ent>              - The specific type of entity
+ *   <EntIter>          - The type of iterator for the entities
+ *   begin              - (I) Iterator to first entity
+ *   end                - (I) Iterator to last entity
+ *   partition          - (I) >  0 - only add elements of this partition to the
+ *                                   zone
+ *                            = -1 - add all elements to the zone (default
+ *                                   value)
+ *
+ *============================================================================*/
+
+template <unsigned DIM>
+template <typename Ent, typename EntIter>
+void MZone<DIM>::add_elements_in_entities
+(EntIter begin, EntIter end, const int partition)
+{
+
+  // Check the Ent and EntIter types
+  // NOTE:  If the compiler sent you here, you're using an invalid entity as a
+  // template parameter and/or and invalid iterator.  See struct
+  // 'ValidEntityIterator' in 'MZone.h' for valid types
+  typedef typename ValidEntityIterator<Ent, typename EntIter::value_type>
+    ::Type Check;
+
+  // Find the neighbours of each vertex, edge, and face
+  for(EntIter itEnt = begin; itEnt != end; ++itEnt) {
+    ParseEntity<Ent>::eval(static_cast<Ent*>(*itEnt), boFaceMap, elemVec,
+                           vertMap, zoneElemConn, partition);
+  }
+
+}
+
+
+/*==============================================================================
+ *
+ * Routine: add_elements_in_entity
+ *
+ * Purpose
+ * =======
+ *
+ *   Adds all (or only those of the right partition) elements in a single entity
+ *   to the zone.
+ *
+ * I/O
+ * ===
+ *
+ *   entity             - (I) Pointer to the entity
+ *   partition          - (I) >  0 - only add elements of this partition to the
+ *                                   zone
+ *                            = -1 - add all elements to the zone (default
+ *                                   value)
+ *
+ * Notes
+ * =====
+ *
+ *   - The entity pointer must be of a specific type (not a base class).
+ *
+ *============================================================================*/
+
+template <unsigned DIM>
+template <typename Ent, typename EntPtr>
+void MZone<DIM>::add_elements_in_entity
+(EntPtr entity, const int partition)
+{
+
+  // Check the type of EntPtr (it must not be GEntity*)
+  // NOTE:  If the compiler sent you here, you're using an invalid entity as a
+  // pointer.  See struct ValidEntityIterator in 'MZone.h' for valid types.
+  typedef typename ValidEntityIterator<Ent, EntPtr>::Type Check;
+
+  // Find the neighbours of each vertex, edge, and face
+  ParseEntity<Ent>::eval(static_cast<Ent*>(entity), boFaceMap, elemVec,
+                         vertMap, zoneElemConn, partition);
+
+}
+
+
+/*==============================================================================
+ *
+ * Routine: zoneData
+ *
+ * Purpose
+ * =======
+ *
+ *   Processes all the elements in a zone.  Stores vertices and element
+ *   connectivity.  Records boundary vertices for use with class
+ *   'MZoneBoundary'.  Both the vertices and elements are ordered with boundary
+ *   vertices/elements first.
+ *
+ *============================================================================*/
+
+template <unsigned DIM>
+int MZone<DIM>::zoneData()
+{
+
+  if(elemVec.size() == 0) return 1;
+
+  // Resize output arrays
+  zoneVertVec.resize(vertMap.size());
+
+//--Label boundary vertices and start building output vector of vertices.
+//--Also record boundary faces that contain a boundary vertex.
+
+  std::vector<MVertex*> faceVertices;
+  unsigned cVert = 0;
+  for(typename BoFaceMap::iterator fMapIt = boFaceMap.begin();
+      fMapIt != boFaceMap.end(); ++fMapIt) {
+    // Get all the vertices on the face
+    FaceTr<FaceT>::getAllFaceVertices
+      (elemVec[fMapIt->second.parentElementIndex].element,
+       fMapIt->second.parentFace, faceVertices);
+    const int nVert = faceVertices.size();
+    for(int iVert = 0; iVert != nVert; ++iVert) {
+      int &index = vertMap[faceVertices[iVert]];
+      // Label the boundary vertices
+      if(index == 0) {
+        zoneVertVec[cVert] = faceVertices[iVert];
+        index = ++cVert;
+      }
+      // Record boundary faces that belong to boundary vertices
+      ZoneVertexData<typename BoFaceMap::const_iterator> &boVertData =
+        boVertMap[faceVertices[iVert]];
+      // const_cast is okay, the keys to 'boFaceMap' will not be changed
+      boVertData.faces.push_back(fMapIt);
+      boVertData.index = index;
+    }
+  }
+  numBoVert = cVert;
+
+//--Label interior vertices and complete output vector of vertices
+
+  const VertexMap::const_iterator vMapEnd = vertMap.end();
+  for(VertexMap::iterator vMapIt = vertMap.begin();
+      vMapIt != vMapEnd; ++vMapIt) {
+    if(vMapIt->second == 0) {  // Vertex in zone interior
+      zoneVertVec[cVert] = vMapIt->first;
+      vMapIt->second = ++cVert;
+    }
+  }
+
+//--Initialize the connectivity array for the various types of elements.  Note
+//--that 'iElemType' is MSH_TYPE-1.
+
+  int numUsedElemType = 0;
+  for(int iElemType = 0; iElemType != MSH_NUM_TYPE; ++iElemType) {
+    if(zoneElemConn[iElemType].numElem > 0) {
+      zoneElemConn[iElemType].connectivity.resize
+        (zoneElemConn[iElemType].numElem*MElement::getInfoMSH(iElemType+1));
+      // Members numBoElem, and iConn should be set to zero via the constructor
+      // or a clear
+    }
+  }
+
+//--The elements are added to the connectivity in two loops.  The first loop
+//--looks for boundary elements.  The second loop adds the remaining interior
+//--elements.
+
+  unsigned cElem = 1;
+  const ElementVec::const_iterator eVecEnd = elemVec.end();
+  for(ElementVec::iterator eVecIt = elemVec.begin(); eVecIt != eVecEnd;
+      ++eVecIt) {
+    // It is sufficient to check the primary vertices to see if this element
+    // is on the boundary
+    const int nPVert = eVecIt->element->getNumPrimaryVertices();
+    for(int iPVert = 0; iPVert != nPVert; ++iPVert) {
+      if(vertMap[eVecIt->element->getVertex(iPVert)] <= numBoVert) {
+        // The element index
+        eVecIt->index = cElem++;
+        // The type of element
+        const int iElemType = eVecIt->element->getTypeForMSH() - 1;
+        // Increment number of boundary elements of this type
+        ++zoneElemConn[iElemType].numBoElem;
+        // Load connectivity for this element type
+        const int nVert = eVecIt->element->getNumVertices();
+        for(int iVert = 0; iVert != nVert; ++iVert) {
+          zoneElemConn[iElemType].add_to_connectivity
+            (vertMap[eVecIt->element->getVertex(iVert)]);
+        }
+        break;
+      }
+    }
+  }
+
+//--Now loop through all elements again and do same thing for all elements with
+//--index set to 0
+
+  for(ElementVec::iterator eVecIt = elemVec.begin(); eVecIt != eVecEnd;
+      ++eVecIt) {
+    if(eVecIt->index == 0) {
+      // The element index
+      eVecIt->index == cElem++;
+      // The type of element
+      const int iElemType = eVecIt->element->getTypeForMSH() - 1;
+      // Load connectivity for this element type
+      const int nVert = eVecIt->element->getNumVertices();
+      for(int iVert = 0; iVert != nVert; ++iVert) {
+        zoneElemConn[iElemType].add_to_connectivity
+          (vertMap[eVecIt->element->getVertex(iVert)]);
+      }
+    }
+  }
+
+//**If we are going to write the boundary element faces, we need to update
+//**.parentElementIndex from "index in elemVec" to "elemVec[iParent].index.
+//**This is the index of the parent element local to the zone.
+
+//--Clean-up for containers that are no longer required.  Only 'boVertMap' is
+//--still required but 'boFaceMap' is also retained because it contains the
+//--faces referenced by 'boVertMap'.
+
+  elemVec.clear();
+  vertMap.clear();
+  return 0;
+
+}
+
+
+/*******************************************************************************
+ *
+ * Template meta-programming classes
+ *
+ ******************************************************************************/
+
+/*--------------------------------------------------------------------*
+ * This class uses traits to determine the element types in an entity
+ * and iterate over them.  A template meta-programming loop is used to
+ * iterate over the types of elements.
+ *--------------------------------------------------------------------*/
+
+//--Entry point
+
+template <typename Ent, int N>
+struct ParseEntity
+{
+  typedef typename EntTr<Ent>::FaceT FaceT;  // The type/dimension of face
+  typedef typename FaceTr<FaceT>::BoFaceMap BoFaceMap;  // The corresponding map
+  typedef typename EntElemTr<Ent, N>::Elem Elem;  // The type of primary element
+  static void eval(const Ent *const entity,
+                   BoFaceMap &boFaceMap,
+                   ElementVec &elemVec,
+                   VertexMap &vertMap,
+                   ElementConnectivity *zoneElemConn,
+                   const int partition)
+  {
+    for(typename EntElemTr<Ent, N>::ConstElementIterator itElem =
+          EntElemTr<Ent, N>::begin(entity);
+        itElem != EntElemTr<Ent, N>::end(entity); ++itElem) {
+      if(partition < 0 || (*itElem)->getPartition() == partition) {
+        // Unique list of elements
+        const unsigned eVecIndex = elemVec.size();
+        elemVec.push_back(ElemData(*itElem));
+        ++zoneElemConn[(*itElem)->getTypeForMSH() - 1].numElem;
+        // Unique list of vertices
+        const int nVert = (*itElem)->getNumVertices();
+        for(int iVert = 0; iVert != nVert; ++iVert)
+          vertMap[(*itElem)->getVertex(iVert)] = 0;  // Unlabelled
+        // Maintain list of (base element) faces with only one bounding element.
+        for(int iFace = 0; iFace != ElemFaceTr<Elem>::numFaceT;
+            ++iFace) {
+          FaceT face = FaceTr<FaceT>::template getFace<Elem>(*itElem, iFace);
+          std::pair<typename BoFaceMap::iterator, bool> insBoFaceMap =
+            boFaceMap.insert(std::pair<FaceT, FaceData>
+                             (face, FaceData(*itElem, iFace, eVecIndex)));
+          if(!insBoFaceMap.second) {
+            // The face already exists and is therefore bounded on both sides
+            // by elements.  Not a boundary face so delete.
+            boFaceMap.erase(insBoFaceMap.first);
+          }
+        }
+      }
+    }
+    // Next element type in the entity
+    ParseEntity<Ent, N-1>::eval(entity, boFaceMap, elemVec, vertMap,
+                                zoneElemConn, partition);
+  }
+};
+
+//--Terminate loop when N = 1
+
+template <typename Ent>
+struct ParseEntity<Ent, 1>
+{
+  typedef typename EntTr<Ent>::FaceT FaceT;  // The type/dimension of face
+  typedef typename FaceTr<FaceT>::BoFaceMap BoFaceMap;  // The corresponding map
+  typedef typename EntElemTr<Ent, 1>::Elem Elem;  // The type of primary element
+  static void eval(const Ent *const entity,
+                   BoFaceMap &boFaceMap,
+                   ElementVec &elemVec,
+                   VertexMap &vertMap,
+                   ElementConnectivity *zoneElemConn,
+                   const int partition)
+  {
+    for(typename EntElemTr<Ent, 1>::ConstElementIterator itElem =
+          EntElemTr<Ent, 1>::begin(entity);
+        itElem != EntElemTr<Ent, 1>::end(entity); ++itElem) {
+      if(partition < 0 || (*itElem)->getPartition() == partition) {
+        // Unique list of elements
+        const unsigned eVecIndex = elemVec.size();
+        elemVec.push_back(ElemData(*itElem));
+        ++zoneElemConn[(*itElem)->getTypeForMSH() - 1].numElem;
+        // Unique list of vertices
+        const int nVert = (*itElem)->getNumVertices();
+        for(int iVert = 0; iVert != nVert; ++iVert)
+          vertMap[(*itElem)->getVertex(iVert)] = 0;  // Unlabelled
+        // Maintain list of (base element) faces with only one bounding element.
+        for(int iFace = 0; iFace != ElemFaceTr<Elem>::numFaceT;
+            ++iFace) {
+          FaceT face = FaceTr<FaceT>::template getFace<Elem>(*itElem, iFace);
+          std::pair<typename BoFaceMap::iterator, bool> insBoFaceMap =
+            boFaceMap.insert(std::pair<FaceT, FaceData>
+                             (face, FaceData(*itElem, iFace, eVecIndex)));
+          if(!insBoFaceMap.second) {
+            // The face already exists and is therefore bounded on both sides
+            // by elements.  Not a boundary face so delete.
+            boFaceMap.erase(insBoFaceMap.first);
+          }
+        }
+      }
+    }
+  }
+};
+
+
+/*******************************************************************************
+ *
+ * Explicit instantiations of class MZone
+ *
+ ******************************************************************************/
+
+//--All the non-template routines in the class
+
+template class MZone<2>;
+template class MZone<3>;
+
+//--Explicit instantiations of the routines for adding elements from a
+//--container of entities
+
+// Vector container
+template void MZone<2>::add_elements_in_entities
+<GFace, std::vector<GEntity*>::const_iterator>
+(std::vector<GEntity*>::const_iterator begin,
+ std::vector<GEntity*>::const_iterator end, const int partition);
+
+template void MZone<3>::add_elements_in_entities
+<GRegion, std::vector<GEntity*>::const_iterator>
+(std::vector<GEntity*>::const_iterator begin,
+ std::vector<GEntity*>::const_iterator end, const int partition);
+
+//--Explicit instantiations of the routines for adding elements from a single
+//--entity
+
+template void MZone<2>::add_elements_in_entity
+<GFace, GFace*>
+(GFace* entity, const int partition);
+
+template void MZone<3>::add_elements_in_entity
+<GRegion, GRegion*>
+(GRegion* entity, const int partition);
diff --git a/Geo/MZone.h b/Geo/MZone.h
new file mode 100644
index 0000000000..4ead3a9f68
--- /dev/null
+++ b/Geo/MZone.h
@@ -0,0 +1,446 @@
+// Gmsh - Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
+//
+// See the LICENSE.txt file for license information. Please report all
+// bugs and problems to <gmsh@geuz.org>.
+//
+// MZone.h - Copyright (C) 2008 S. Guzik, C. Geuzaine, J.-F. Remacle
+#ifndef _MZONE_H_
+#define _MZONE_H_
+
+
+/*******************************************************************************
+ *
+ * The classes in this file are used to define and generate representations of
+ * zones.
+ *
+ ******************************************************************************/
+
+#include <map>
+#include <vector>
+#include "MElement.h"
+#include "GFace.h"
+#include "GRegion.h"
+#include "MEdgeHash.h"
+#include "MFaceHash.h"
+#include "CustomContainer.h"
+
+// #define HAVE_HASH_MAP
+
+#if defined(HAVE_HASH_MAP)
+#include "HashMap.h"
+#endif
+
+
+/*==============================================================================
+ * Forward declarations
+ *============================================================================*/
+
+template <unsigned DIM>
+class MZoneBoundary;
+
+
+/*==============================================================================
+ * Required types
+ *============================================================================*/
+
+//--Record of unique elements
+
+struct ElemData
+{
+  MElement *element;
+  int index;
+  ElemData(MElement *const _element)
+    : element(_element), index(0)
+  { }
+};
+
+typedef std::vector<ElemData> ElementVec;
+
+//--Record of unique vertices
+
+typedef std::map<MVertex*, int, std::less<MVertex*> > VertexMap;
+
+//--Data for each face.  Ultimately, only faces on the boundary of the zone are
+//--stored.  Value type for 'BoFaceMap'.
+
+struct FaceData
+{
+  MElement *parentElement;
+  int parentFace;
+  int parentElementIndex;
+  int faceIndex;
+  FaceData(MElement *const _parentElement, const int _parentFace,
+           const int _parentElementIndex)
+    : parentElement(_parentElement), parentFace(_parentFace),
+    parentElementIndex(_parentElementIndex), faceIndex(0)
+  { }
+};
+
+//--Provides information and boundary faces connected to a vertex.  Value type
+//--for 'BoVertexMap'
+
+template <typename BFMapIt>
+struct ZoneVertexData
+{
+  CCon::FaceVector<BFMapIt> faces; // Vector optimized for storing faces
+  int index;
+};
+
+/*--------------------------------------------------------------------*
+ * User I/O
+ *--------------------------------------------------------------------*/
+
+struct ElementConnectivity
+{
+  std::vector<int> connectivity;
+  int numElem;
+  int numBoElem;
+  int iConn;
+  // Constructor
+  ElementConnectivity()
+    : numElem(0), numBoElem(0), iConn(0)
+  { }
+  // Member functions
+  void add_to_connectivity(const int i)
+  {
+    connectivity[iConn++] = i;
+  }
+  void clear()
+  {
+    connectivity.clear();
+    numElem = 0;
+    numBoElem = 0;
+    iConn = 0;
+  }
+};
+
+//--Output type for zone element connectivity
+
+typedef std::vector<ElementConnectivity> ElementConnectivityVec;
+
+//--Output type for vertices in the zone
+
+typedef std::vector<MVertex*> VertexVec;
+
+
+/*==============================================================================
+ * Traits classes - that just perform compile-time checks for valid argument
+ * types
+ *============================================================================*/
+
+//--This traits class checks for valid types of entities and iterators.
+
+// NOTE:  If the compiler sent you here, you're using an invalid entity and/or
+// invalid iterator.  Valid pairs are shown below.
+template <typename Ent, typename Iter> struct ValidEntityIterator;
+template <> struct ValidEntityIterator<GFace, GFace*>
+{ typedef void Type; };
+template <> struct ValidEntityIterator<GFace, GEntity*>
+{ typedef void Type; };
+template <> struct ValidEntityIterator<GRegion, GRegion*>
+{ typedef void Type; };
+template <> struct ValidEntityIterator<GRegion, GEntity*>
+{ typedef void Type; };
+
+
+/*==============================================================================
+ * Traits classes - that return information about a type
+ *============================================================================*/
+
+//--Traits based on the dimension
+
+template <unsigned DIM> struct DimTr;
+template <> struct DimTr<2>
+{
+  typedef MEdge FaceT;
+  typedef GFace EntityT;
+};
+template <> struct DimTr<3>
+{
+  typedef MFace FaceT;
+  typedef GRegion EntityT;
+};
+
+//--This traits class describes the number of dimension-based 'FaceT' in an
+//--primary element type
+
+template <typename Elem> struct ElemFaceTr;
+template <> struct ElemFaceTr<MTriangle>    { enum { numFaceT =  3 }; };
+template <> struct ElemFaceTr<MQuadrangle>  { enum { numFaceT =  4 }; };
+template <> struct ElemFaceTr<MTetrahedron> { enum { numFaceT =  4 }; };
+template <> struct ElemFaceTr<MHexahedron>  { enum { numFaceT =  6 }; };
+template <> struct ElemFaceTr<MPrism>       { enum { numFaceT =  5 }; };
+template <> struct ElemFaceTr<MPyramid>     { enum { numFaceT =  5 }; };
+
+//--This traits class gives the number of element types in entity Ent
+
+template <typename Ent> struct EntTr;
+template <> struct EntTr<GFace>
+{
+  typedef MEdge FaceT;
+  enum { numElemTypes = 2 };
+};
+template <> struct EntTr<GRegion>
+{
+  typedef MFace FaceT;
+  enum { numElemTypes = 4 };
+};
+
+//--This traits class gives iterator types and begin and end iterators for
+//--element type number N in entity Ent.
+
+template <typename Ent, int N> struct EntElemTr;
+template <> struct EntElemTr<GFace, 1> {
+  typedef MQuadrangle Elem;
+  typedef std::vector<MQuadrangle*>::const_iterator ConstElementIterator;
+  static ConstElementIterator begin(const GFace *const gFace)
+  {
+    return gFace->quadrangles.begin();
+  }
+  static ConstElementIterator end(const GFace *const gFace)
+  {
+    return gFace->quadrangles.end();
+  }
+};
+template <> struct EntElemTr<GFace, 2> {
+  typedef MTriangle Elem;
+  typedef std::vector<MTriangle*>::const_iterator ConstElementIterator;
+  static ConstElementIterator begin(const GFace *const gFace)
+  {
+    return gFace->triangles.begin();
+  }
+  static ConstElementIterator end(const GFace *const gFace)
+  {
+    return gFace->triangles.end();
+  }
+};
+template <> struct EntElemTr<GRegion, 1> {
+  typedef MPyramid Elem;
+  typedef std::vector<MPyramid*>::const_iterator ConstElementIterator;
+  static ConstElementIterator begin(const GRegion *const gRegion)
+  {
+    return gRegion->pyramids.begin();
+  }
+  static ConstElementIterator end(const GRegion *const gRegion)
+  {
+    return gRegion->pyramids.end();
+  }
+};
+template <> struct EntElemTr<GRegion, 2> {
+  typedef MPrism Elem;
+  typedef std::vector<MPrism*>::const_iterator ConstElementIterator;
+  static ConstElementIterator begin(const GRegion *const gRegion)
+  {
+    return gRegion->prisms.begin();
+  }
+  static ConstElementIterator end(const GRegion *const gRegion)
+  {
+    return gRegion->prisms.end();
+  }
+};
+template <> struct EntElemTr<GRegion, 3> {
+  typedef MHexahedron Elem;
+  typedef std::vector<MHexahedron*>::const_iterator ConstElementIterator;
+  static ConstElementIterator begin(const GRegion *const gRegion)
+  {
+    return gRegion->hexahedra.begin();
+  }
+  static ConstElementIterator end(const GRegion *const gRegion)
+  {
+    return gRegion->hexahedra.end();
+  }
+};
+template <> struct EntElemTr<GRegion, 4> {
+  typedef MTetrahedron Elem;
+  typedef std::vector<MTetrahedron*>::const_iterator ConstElementIterator;
+  static ConstElementIterator begin(const GRegion *const gRegion)
+  {
+    return gRegion->tetrahedra.begin();
+  }
+  static ConstElementIterator end(const GRegion *const gRegion)
+  {
+    return gRegion->tetrahedra.end();
+  }
+};
+
+//--Traits/policy class based on dimension of a face.
+
+template <typename FaceT> struct FaceTr;
+template <> struct FaceTr<MEdge> {
+  typedef std::map<MEdge, FaceData, Less_Edge> BoFaceMap;
+  template <typename Elem>
+  static MEdge getFace(Elem *const element, const int iFace)
+  {
+    return element->getEdge(iFace);  // Primary element
+  }
+  static void getAllFaceVertices(MElement *const element, const int iFace,
+                                 std::vector<MVertex*> &v)
+  {
+    element->getEdgeVertices(iFace, v);  // Any element
+  }
+};
+template <> struct FaceTr<MFace> {
+  typedef std::map<MFace, FaceData, Less_Face> BoFaceMap;
+  template <typename Elem>
+  static MFace getFace(Elem *const element, const int iFace)
+  {
+    return element->getFace(iFace);  // Primary element
+  }
+  static void getAllFaceVertices(MElement *const element, const int iFace,
+                                 std::vector<MVertex*> &v)
+  {
+    element->getFaceVertices(iFace, v);  // Any element
+  }
+};
+
+
+/*******************************************************************************
+ *
+ * class: MZone
+ *
+ * Purpose
+ * =======
+ *
+ *   Generates a definition of a zone based on entities and/or partitions.
+ *
+ *   Template parameters:
+ *     DIM              - dimension of the problem
+ *
+ * Notes
+ * =====
+ *
+ *   - explicitly instantiated in 'MZone.cpp'
+ *
+ ******************************************************************************/
+
+template <unsigned DIM>
+class MZone
+{
+
+
+/*==============================================================================
+ * Internal types
+ *============================================================================*/
+
+ private:
+  typedef typename DimTr<DIM>::FaceT FaceT;
+  typedef typename FaceTr<FaceT>::BoFaceMap BoFaceMap;
+  typedef typename std::map<const MVertex*,
+    ZoneVertexData<typename BoFaceMap::const_iterator>,
+    std::less<const MVertex*> > BoVertexMap;
+
+
+/*==============================================================================
+ * Member functions
+ *============================================================================*/
+
+ public:
+
+//--Default constructor.
+
+  MZone()
+    : numBoVert(0)
+  {
+    elemVec.reserve(8192);
+  }
+
+/*--------------------------------------------------------------------*
+ * Elements added from entities.
+ * Note: It is much easier to keep these in the .cpp file but that
+ *   requries explicit instantiations for each Ent and EntIter.
+ *   Currently, instantiations only exist for containers of type:
+ *     vector
+ *   More can be added as required at the end of MZone.cpp
+ *--------------------------------------------------------------------*/
+
+//--Add all elements in a container of entities.  The specific type of entity
+//--is not known and must be specified as parameter 'Ent'.
+
+  template <typename Ent, typename EntIter>
+  void add_elements_in_entities(EntIter begin, EntIter end,
+                                const int partition = -1);
+
+//--Add elements in a single entity.
+
+  template <typename Ent, typename EntPtr>
+  void add_elements_in_entity(EntPtr entity, const int partition = -1);
+  
+
+/*--------------------------------------------------------------------*
+ * Reset the database
+ *--------------------------------------------------------------------*/
+
+  void clear() 
+  {
+    elemVec.clear();
+    vertMap.clear();
+    boFaceMap.clear();
+    boVertMap.clear();
+    zoneVertVec.clear();
+    for(int i = 0; i != MSH_NUM_TYPE; ++i) {
+      zoneElemConn[i].clear();
+    }
+    numBoVert = 0;
+  }
+
+/*--------------------------------------------------------------------*
+ * Process/query the zone - only after all elements have been added!
+ *--------------------------------------------------------------------*/
+
+//--Compute the zone data
+
+  int zoneData();
+
+//--Total number of elements
+
+  int totalElements() const 
+  {
+    int numElem = 0;
+    for(int iElemType = 0; iElemType != MSH_NUM_TYPE; ++iElemType)
+      numElem += zoneElemConn[iElemType].numElem;
+    return numElem;
+  }
+
+//--Number of element types
+
+  int numElementTypes() const
+  {
+    int numElemType = 0;
+    for(int iElemType = 0; iElemType != MSH_NUM_TYPE; ++iElemType)
+      if(zoneElemConn[iElemType].numElem > 0) ++numElemType;
+    return numElemType;
+  }
+
+
+/*==============================================================================
+ * Member data
+ *============================================================================*/
+
+ private:
+
+//--Data members
+
+  ElementVec elemVec;                   // Set of unique elements
+  VertexMap vertMap;                    // Set of unique vertices and associated
+                                        // numbers in the zone
+  BoFaceMap boFaceMap;                  // Map of boundary faces
+  BoVertexMap boVertMap;                // Map of boundary vertices
+
+ public:
+
+  // I/O
+  VertexVec zoneVertVec;
+  ElementConnectivity zoneElemConn[MSH_NUM_TYPE];
+                                        // Connectivity for each type of element
+  int numBoVert;
+
+
+/*==============================================================================
+ * Friends
+ *============================================================================*/
+
+  friend class MZoneBoundary<DIM>;
+
+};
+
+#endif
diff --git a/Geo/MZoneBoundary.cpp b/Geo/MZoneBoundary.cpp
new file mode 100644
index 0000000000..7027f31917
--- /dev/null
+++ b/Geo/MZoneBoundary.cpp
@@ -0,0 +1,902 @@
+// Gmsh - Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
+//
+// See the LICENSE.txt file for license information. Please report all
+// bugs and problems to <gmsh@geuz.org>.
+//
+// MZoneBoundary.cpp - Copyright (C) 2008 S. Guzik, C. Geuzaine, J.-F. Remacle
+
+#include <iostream> // DBG
+#include <limits> // ?
+
+#include "MZoneBoundary.h"
+#include "Message.h"
+
+
+/*******************************************************************************
+ *
+ * class BCPatchIndex
+ *
+ * Purpose
+ * =======
+ *
+ *   Keeps track of the BC patch index for each entity index.  Many entities
+ *   may share the same patch index and this class enables multiple entities to
+ *   share the same index.
+ *
+ ******************************************************************************/
+
+class BCPatchIndex
+{
+  struct PatchData 
+  {
+    int index;
+    std::vector<int> eTagVec;
+    PatchData(int eTag1)
+      : index(eTag1), eTagVec(1, eTag1)
+    { }
+  };
+
+  typedef std::list<PatchData> PatchDataList;
+  typedef PatchDataList::iterator PatchDataListIt;
+
+  typedef std::map<int, PatchDataListIt> PatchMap;
+  typedef PatchMap::value_type PatchMapVal;
+  typedef PatchMap::iterator PatchMapIt;
+  typedef std::pair<PatchMapIt, bool> PatchMapIns;
+
+  // Data
+  PatchDataList patchData;
+  PatchMap patch;
+  bool sharedPatch;
+
+  PatchMapIt _add(const int eTag)
+  {
+    PatchMapIns insPatch = patch.insert(PatchMapVal(eTag, PatchDataListIt()));
+    if(insPatch.second) {
+      insPatch.first->second =
+        patchData.insert(patchData.end(), PatchData(eTag));
+    }
+    return insPatch.first;
+  }
+
+ public:
+
+  BCPatchIndex()
+    : sharedPatch(false)
+  { }
+
+  // Add an entity tag
+  void add(const int eTag) { _add(eTag); }
+
+  // Add two entity tags which must be given the same patch index
+  void addPair(const int eTag1, const int eTag2)
+  {
+    sharedPatch = true;
+    PatchMapIt patch1 = _add(eTag1);
+    PatchMapIt patch2 = _add(eTag2);
+    if(patch1->second != patch2->second) {
+      PatchData &PD1 = *(patch1->second);
+      PatchData &PD2 = *(patch2->second);
+      const int nTag = PD2.eTagVec.size();
+      for(int iTag = 0; iTag != nTag; ++iTag) {
+        // Move tag from PD2 to PD1
+        const int tag = PD2.eTagVec[iTag];
+        PD1.eTagVec.push_back(tag);
+        // Update value in 'patch' for this tag
+        if(tag != eTag2) patch[tag] = patch1->second;
+      }
+      patchData.erase(patch2->second);
+      patch2->second = patch1->second;
+    }
+  }
+
+  // Once all entity tags have been added, generate patch indices
+  void generatePatchIndices()
+  {
+    if(sharedPatch) {  // Don't renumber if not shared to preserve entity
+                       // numbers
+      int c = 1;
+      for(PatchDataListIt pDIt = patchData.begin(); pDIt != patchData.end();
+          ++pDIt) pDIt->index = c++;
+    }
+  }
+
+  // Get the patch index for an entity (generate patch indices first)
+  int getIndex(const int eTag) { return patch[eTag]->index; }
+
+};
+
+
+/*******************************************************************************
+ *
+ * Routine updateBoVec
+ *
+ * Purpose
+ * =======
+ *
+ *   This routine basically tries to determine the normal for each boundary
+ *   vertex.  Attempts are made to get the normal from the geometry.  If this
+ *   fails, an average of the relavent mesh elements is used.
+ *
+ *   Specializations exist for 2-D and 3-D.
+ *
+ *   If a vertex lies on a lower dimension entity (GVertex in 2-D, GVertex or
+ *   GEdge in 3-D), multiple BC patches will be written in case the normals are
+ *   different for the multiple GEdges (2-D) or GFaces (3-D) that connect to the
+ *   vertex.  In 2-D, if the angle between the normals is less than (1deg?), the
+ *   patches are merged.  This has not yet been implemented in 3D (see notes in
+ *   the 3D specialization below).
+ *
+ *   Normally, BC patches are numbered based on the entities.  If entities are
+ *   merged, class 'BCPatchIndex' is used to give the same patch index to
+ *   multiple entities.
+ *
+ * I/O
+ * ===
+ *
+ *   vertex             - (I) vertex to add to zoneBoVec and to find the normal
+ *                            at
+ *   zoneIndex          - (I)
+ *   vertIndex          - (I) index of the vertex in the zone
+ *   faces              - (I) all boundary faces connected to the vertex
+ *   zoneBoVec          - (O) updated with domain boundary information for the
+ *                            vertex
+ *   patch              - (O) record BC patch index for an entity
+ *   warnNormFromElem   - (I) wether a warning about obtaining normals from
+ *                            elements has already been given. 
+ *                        (O) set to true if warning given in this call
+ *
+ * Notes
+ * =====
+ *
+ *   - 'getNormalFromElements' is at the bottom.  In that case, normals are
+ *     determined from the elements and not the geometry.
+ *   - 'vertex' is the key to 'globalBoVertMap'.  It will not change so the
+ *     const_cast used below is okay.
+ *
+ ******************************************************************************/
+
+template <unsigned DIM, typename FaceT>
+void updateBoVec
+(const MVertex *const vertex, const int zoneIndex, const int vertIndex,
+ const CCon::FaceVector<typename MZoneBoundary<DIM>::template
+ GlobalVertexData<FaceT>::FaceDataB> &faces, ZoneBoVec &zoneBoVec,
+ BCPatchIndex &patch, bool &warnNormFromElem);
+
+
+/*******************************************************************************
+ *
+ * Routine updateBoVec - SPECIALIZED for 2D
+ *
+ ******************************************************************************/
+
+//--"Contained" routines
+
+
+/*==============================================================================
+ *
+ * Routine edge_normal
+ *
+ * Purpose
+ * =======
+ *
+ *   Gives the edge normal.  Use only by routine updateBoVec<2, MEdge>.
+ *
+ * I/O
+ * ===
+ *
+ *   vertex             - (I) The vertex to find the normal at
+ *   edge               - (I) The geometry edge upon which the vertex resides
+ *                            and from which the normal will be determined
+ *   faces              - (I) All faces on the boundary connected to 'vertex'
+ *   boNormal           - (O) The normal to the boundary face (edge in 2D)
+ *   returns            - (O) 0 - success
+ *                            1 - parFromPoint() failed
+ *
+ *============================================================================*/
+
+int edge_normal
+(const MVertex *const vertex, const GEdge *const gEdge,
+ const CCon::FaceVector<MZoneBoundary<2>::GlobalVertexData<MEdge>::FaceDataB>
+ &faces, SVector3 &boNormal)
+{
+  const double par = gEdge->parFromPoint(vertex->point());
+  if(par == std::numeric_limits<double>::max()) return 1;
+
+  const SVector3 tangent(gEdge->firstDer(par));
+                                        // Tangent to the boundary face
+  SPoint3 interior(0., 0., 0.);         // An interior point
+  SVector3 meshPlaneNormal(0.);         // This normal is perpendicular to the
+                                        // plane of the mesh
+
+  // The interior point and mesh plane normal are computed from all elements.
+  const int nFace = faces.size();
+  for(int iFace = 0; iFace != nFace; ++iFace) {
+    interior += faces[iFace].parentElement->barycenter();
+    // Make sure all the planes go in the same direction
+    //**Required?
+    SVector3 mpnt = faces[iFace].parentElement->getFace(0).normal();
+    if(dot(mpnt, meshPlaneNormal) < 0.) mpnt.negate();
+    meshPlaneNormal += mpnt;
+  }
+  interior /= nFace;
+  // Normal to the boundary edge (but unknown direction)
+  boNormal = crossprod(tangent, meshPlaneNormal);
+  boNormal.normalize();
+  // Direction vector from vertex to interior (inwards).  The normal should
+  // point in the same direction.
+  if(dot(boNormal, SVector3(vertex->point(), interior)) < 0.)
+    boNormal.negate();
+  return 0;
+}
+
+
+/*==============================================================================
+ *
+ * Routine updateBoVec<2, MEdge>
+ *
+ *============================================================================*/
+
+template <>
+void updateBoVec<2, MEdge>
+(const MVertex *const vertex, const int zoneIndex, const int vertIndex,
+ const CCon::FaceVector<MZoneBoundary<2>::GlobalVertexData<MEdge>::FaceDataB>
+ &faces, ZoneBoVec &zoneBoVec, BCPatchIndex &patch, bool &warnNormFromElem)
+{
+  SVector3 boNormal;                    // Normal to the boundary face (edge in
+                                        // 2D)
+
+  GEntity *const ent = vertex->onWhat();
+  if(ent == 0) {
+    goto getNormalFromElements;
+    // No entity: try to find a normal from the faces
+  }
+  else {
+    switch(ent->dim()) {
+    case 0:
+
+/*--------------------------------------------------------------------*
+ * In this case, there are possibly two GEdges from this zone
+ * connected to the vertex.  If the angle between the two edges is
+ * significant, the BC patch will be split and this vertex will be
+ * written in both patches with different normals.  If the angle is
+ * small, or only one GEdge belongs to this zone, then the vertex will
+ * only be written to one patch.
+ *--------------------------------------------------------------------*/
+
+      {
+        // Get the edge entities that are connected to the vertex
+        std::list<GEdge*> gEdgeList = ent->edges();
+        // Find edge entities that are connected to elements in the zone
+        std::vector<GEdge*> useGEdge;
+        const int nFace = faces.size();
+        for(int iFace = 0; iFace != nFace; ++iFace) {
+          if(zoneIndex == faces[iFace].zoneIndex) {
+            MEdge mEdge = faces[iFace].parentElement->getEdge
+              (faces[iFace].parentFace);
+            // Get the other vertex on the mesh edge.
+            MVertex *vertex2 = mEdge.getMinVertex();
+            if(vertex2 == vertex) vertex2 = mEdge.getMaxVertex();
+            // Check if the entity associated with vertex2 is a line and
+            // is also connected to vertex.  If so, add it to the container of
+            // edge entities that will be used to determine the normal
+            GEntity *const ent2 = vertex2->onWhat();
+            if(ent2->dim() == 1) useGEdge.push_back(static_cast<GEdge*>(ent2));
+          }
+        }
+
+//--'useGEdge' now contains the edge entities that will be used to determine
+//--the normals
+
+        switch(useGEdge.size()) {
+
+        case 0:
+          goto getNormalFromElements;
+
+        case 1:
+          {
+            const GEdge *const gEdge =
+              static_cast<const GEdge*>(useGEdge[0]);
+            SVector3 boNormal;
+            if(edge_normal(vertex, gEdge, faces, boNormal))
+               goto getNormalFromElements;
+            zoneBoVec.push_back(VertexBoundary(zoneIndex, gEdge->tag(),
+                                               boNormal,
+                                               const_cast<MVertex*>(vertex),
+                                               vertIndex));
+            patch.add(gEdge->tag());
+          }
+          break;
+
+        case 2:
+          {
+            // Get the first normal
+            const GEdge *const gEdge1 =
+              static_cast<const GEdge*>(useGEdge[0]);
+            SVector3 boNormal1;
+            if(edge_normal(vertex, gEdge1, faces, boNormal1))
+              goto getNormalFromElements;
+
+            // Get the second normal
+            const GEdge *const gEdge2 =
+              static_cast<const GEdge*>(useGEdge[1]);
+            SVector3 boNormal2;
+            if(edge_normal(vertex, gEdge2, faces, boNormal2))
+              goto getNormalFromElements;
+
+            if(dot(boNormal1, boNormal2) < 0.98) {
+
+//--Write 2 separate patches
+
+              zoneBoVec.push_back(VertexBoundary(zoneIndex, gEdge1->tag(),
+                                                 boNormal1,
+                                                 const_cast<MVertex*>(vertex),
+                                                 vertIndex));
+              patch.add(gEdge1->tag());
+              zoneBoVec.push_back(VertexBoundary(zoneIndex, gEdge2->tag(),
+                                                 boNormal2,
+                                                 const_cast<MVertex*>(vertex),
+                                                 vertIndex));
+              patch.add(gEdge2->tag());
+            }
+            else {
+
+//--Write one patch
+
+              boNormal1 += boNormal2;
+              boNormal1 *= 0.5;
+              zoneBoVec.push_back(VertexBoundary(zoneIndex, gEdge1->tag(),
+                                                 boNormal1,
+                                                 const_cast<MVertex*>(vertex),
+                                                 vertIndex));
+              patch.addPair(gEdge1->tag(), gEdge2->tag());
+            }
+          }
+          break;
+
+        default:
+          Msg::Error("Internal error based on 2-D boundary assumptions (file "
+                     "\'MZoneBoundary.cpp').  Boundary vertices may be "
+                     "incorrect");
+          break;
+        }
+      }
+      break;
+    case 1:
+
+/*--------------------------------------------------------------------*
+ * The vertex exists on an edge and belongs to only 1 BC patch.
+ *--------------------------------------------------------------------*/
+
+      {
+        SVector3 boNormal;
+        if(edge_normal(vertex, static_cast<const GEdge*>(ent), faces, boNormal))
+          goto getNormalFromElements;
+        zoneBoVec.push_back(VertexBoundary(zoneIndex, ent->tag(), boNormal,
+                                           const_cast<MVertex*>(vertex),
+                                           vertIndex));
+        patch.add(ent->tag());
+      }
+      break;
+
+    default:
+      goto getNormalFromElements;
+    }
+  }
+  return;
+
+ getNormalFromElements: ;
+
+/*--------------------------------------------------------------------*
+ * Geometry information cannot be used - generate normals from the
+ * elements
+ *--------------------------------------------------------------------*/
+
+  {
+    if(warnNormFromElem) {
+      Msg::Warning("Some or all boundary normals were determined from mesh "
+                   "elements instead of from the geometry");
+      warnNormFromElem = false;
+    }
+    // The mesh plane normal is computed from all elements attached to the
+    // vertex
+    SVector3 meshPlaneNormal(0.);       // This normal is perpendicular to the
+                                        // plane of the mesh
+    const int nFace = faces.size();
+    for(int iFace = 0; iFace != nFace; ++iFace) {
+      // Make sure all the planes go in the same direction
+      //**Required?
+      SVector3 mpnt = faces[iFace].parentElement->getFace(0).normal();
+      if(dot(mpnt, meshPlaneNormal) < 0.) mpnt.negate();
+      meshPlaneNormal += mpnt;
+    }
+    // Sum the normals from each element.  The tangent is computed from all
+    // faces in the zone attached to the vertex and is weighted by the length of
+    // the edge.  Each tangent has to be converted independently into an
+    // inwards-pointing normal.
+    SVector3 boNormal(0.);
+    for(int iFace = 0; iFace != nFace; ++iFace) {
+      if(zoneIndex == faces[iFace].zoneIndex) {
+        const SVector3 tangent = faces[iFace].parentElement->getEdge
+          (faces[iFace].parentFace).tangent();
+        // Normal to the boundary (unknown direction)
+        SVector3 bnt = crossprod(tangent, meshPlaneNormal);
+        // Inwards normal
+        const SVector3 inwards(vertex->point(),
+                               faces[iFace].parentElement->barycenter());
+        if(dot(bnt, inwards) < 0.) bnt.negate();
+        boNormal += bnt;
+      }
+    }
+    boNormal.normalize();
+    zoneBoVec.push_back(VertexBoundary(zoneIndex, 0, boNormal,
+                                       const_cast<MVertex*>(vertex),
+                                       vertIndex));
+    patch.add(0);
+  }
+}
+
+
+/*******************************************************************************
+ *
+ * Routine updateBoVec - SPECIALIZED for 3D
+ *
+ * Notes
+ * =====
+ *
+ *   - Merging two entities into a single patch is more difficult in 3D because,
+ *     while in 2D, edge entities are separated by a single vertex, in 3D
+ *     entities are separated by a number of vertices, each of which may force
+ *     splitting of the patch.  Merging of entities into a single patch has not
+ *     yet been implemented, but the best way to accomplish this is probably to
+ *     split all entities into separate patches and then record which pairs of
+ *     entities cannot be merged.  Later, all others and be merged and
+ *     'zoneBoVec' updated accordingly.
+ *
+ ******************************************************************************/
+
+template <>
+void updateBoVec<3, MFace>
+(const MVertex *const vertex, const int zoneIndex, const int vertIndex,
+ const CCon::FaceVector<MZoneBoundary<3>::GlobalVertexData<MFace>::FaceDataB>
+ &faces, ZoneBoVec &zoneBoVec, BCPatchIndex &patch, bool &warnNormFromElem)
+{
+  GEntity *const ent = vertex->onWhat();
+  if(ent == 0) {
+    goto getNormalFromElements;
+    // No entity: try to find a normal from the faces
+  }
+  else {
+    switch(ent->dim()) {
+    case 0:
+    case 1:
+
+/*--------------------------------------------------------------------*
+ * In this case, there are possibly multiple GFaces from this zone
+ * connected to the vertex.  One patch for each GFace will be written.
+ *--------------------------------------------------------------------*/
+
+      {
+
+//--Get a list of face entities connected to 'ent'
+
+        std::list<GFace*> gFaceList;
+        switch(ent->dim()) {
+        case 0:
+          {
+            std::list<GEdge*> gEdgeList = ent->edges();
+            std::list<GFace*> gFaceList;
+            for(std::list<GEdge*>::const_iterator gEIt = gEdgeList.begin();
+                gEIt != gEdgeList.end(); ++gEIt) {
+              std::list<GFace*> alist = (*gEIt)->faces();
+              gFaceList.splice(gFaceList.end(), alist);
+              // gcc-4.3 won't accept the following
+//               gFaceList.splice(gFaceList.end(), (*gEIt)->faces());
+            }
+            // Remove duplicates
+            gFaceList.sort();
+            gFaceList.unique();
+          }
+          break;
+        case 1:
+          gFaceList = ent->faces();
+          break;
+        }
+
+//--Get a list of face entities connected to the 'vertex' that are also in the
+//--zone
+
+        std::list<const GFace*> useGFace;
+        std::list<GEdge*> gEdgeList;
+        const int nFace = faces.size();
+        for(int iFace = 0; iFace != nFace; ++iFace) {
+          if(zoneIndex == faces[iFace].zoneIndex) {
+            bool matchedFace = false;
+            MFace mFace = faces[iFace].parentElement->getFace
+              (faces[iFace].parentFace);
+            const int nVOnF = mFace.getNumVertices();
+            int vertexOnF;              // The index of 'vertex' in the face
+            for(int iVOnF = 0; iVOnF != nVOnF; ++iVOnF) {
+              const MVertex *const vertex2 = mFace.getVertex(iVOnF);
+              if(vertex == vertex2) vertexOnF = iVOnF;
+              else {
+                const GEntity *const ent2 = vertex2->onWhat();
+                if(ent2->dim() == 2) {
+                  matchedFace = true;
+                  useGFace.push_back(static_cast<const GFace*>(ent2));
+                  break;
+                }
+              }
+            }
+            // Triangle MElements are a total P.I.T.A.:
+            // - If the original 'ent' is a vertex, one MVertex can be on the
+            //   GVertex, and the other two on GEdges, and then the MElement is
+            //   still on the GFace.
+            // - If the original 'ent' is an edge, one MVertex can be on the
+            //   original GEdge, another on a GVertex, and the final on another
+            //   GEdge, and then the MElement is still on the GFace.  There is
+            //   also the unlikely case where the two other MVertex are both on
+            //   edges ... and the MElement is still on the GFace.
+            if(!matchedFace && nVOnF == 3) {
+              const MVertex *vertex2;
+              const MVertex *vertex3;
+              switch(vertexOnF) {
+              case 0:
+                vertex2 = mFace.getVertex(1);
+                vertex3 = mFace.getVertex(2);
+                break;
+              case 1:
+                vertex2 = mFace.getVertex(0);
+                vertex3 = mFace.getVertex(2);
+                break;
+              case 2:
+                vertex2 = mFace.getVertex(0);
+                vertex3 = mFace.getVertex(1);
+                break;
+              }
+              const GEntity *const ent2 = vertex2->onWhat();
+              const GEntity *const ent3 = vertex3->onWhat();
+              if(ent2->dim() == 1 && ent3->dim() == 1) {
+                // Which GFace is bounded by edges ent2 and ent3?
+                for(std::list<GFace*>::const_iterator gFIt = gFaceList.begin();
+                    gFIt != gFaceList.end(); ++gFIt) {
+                  gEdgeList = (*gFIt)->edges();
+                  if((std::find(gEdgeList.begin(), gEdgeList.end(), ent2)
+                      != gEdgeList.end()) &&
+                     (std::find(gEdgeList.begin(), gEdgeList.end(), ent3)
+                      != gEdgeList.end())) {
+                    // Edges ent2 and ent3 bound this face
+                    useGFace.push_back(*gFIt);
+                    break;
+                  }
+                }
+              }
+              else if(ent->dim() == 1 && (ent2->dim() + ent3->dim() == 1)) {
+                const GEntity *entCmp;
+                if(ent2->dim() == 1) entCmp = ent2;
+                else entCmp = ent3;
+                // Which GFace is bounded by entCmp
+                for(std::list<GFace*>::const_iterator gFIt = gFaceList.begin();
+                    gFIt != gFaceList.end(); ++gFIt) {
+                  gEdgeList = (*gFIt)->edges();
+                  if(std::find(gEdgeList.begin(), gEdgeList.end(), entCmp)
+                     != gEdgeList.end()) {
+                    // Edge entCmp and the original edge bound this face
+                    useGFace.push_back(*gFIt);
+                    break;
+                  }
+                }
+              }
+            }  // Stupid triangles
+          }  // End if face in zone
+        }  // End loop over faces
+        // Duplicates are a possibility, remove
+        useGFace.sort();
+        useGFace.unique();
+
+//--'useGFace' now contains the face entities that connect to vertex.  A BC
+//--patch will be written for each of them.
+
+        for(std::list<const GFace*>::const_iterator gFIt = useGFace.begin();
+            gFIt != useGFace.end(); ++gFIt) {
+          const SPoint2 par = (*gFIt)->parFromPoint(vertex->point());
+          if(par.x() == std::numeric_limits<double>::max())
+            goto getNormalFromElements;  // :P  After all that!
+
+          SVector3 boNormal = (*gFIt)->normal(par);
+          SPoint3 interior(0., 0., 0.);
+          const int nFace = faces.size();
+          for(int iFace = 0; iFace != nFace; ++iFace)
+            interior += faces[iFace].parentElement->barycenter();
+          interior /= nFace;
+          if(dot(boNormal, SVector3(vertex->point(), interior)) < 0.)
+            boNormal.negate();
+
+          zoneBoVec.push_back(VertexBoundary(zoneIndex, (*gFIt)->tag(),
+                                             boNormal,
+                                             const_cast<MVertex*>(vertex),
+                                             vertIndex));
+          patch.add((*gFIt)->tag());
+        }
+      }
+      break;
+    case 2:
+
+/*--------------------------------------------------------------------*
+ * The vertex exists on a face and belongs to only 1 BC patch.
+ *--------------------------------------------------------------------*/
+
+      {
+        const GFace *const gFace = static_cast<const GFace*>(ent);
+        const SPoint2 par = gFace->parFromPoint(vertex->point());
+        if(par.x() == std::numeric_limits<double>::max())
+          goto getNormalFromElements;
+
+        SVector3 boNormal = static_cast<const GFace*>(ent)->normal(par);
+        SPoint3 interior(0., 0., 0.);
+        const int nFace = faces.size();
+        for(int iFace = 0; iFace != nFace; ++iFace)
+          interior += faces[iFace].parentElement->barycenter();
+        interior /= nFace;
+        if(dot(boNormal, SVector3(vertex->point(), interior)) < 0.)
+          boNormal.negate();
+
+        zoneBoVec.push_back(VertexBoundary(zoneIndex, gFace->tag(), boNormal,
+                                           const_cast<MVertex*>(vertex),
+                                           vertIndex));
+        patch.add(gFace->tag());
+      }
+      break;
+    default:
+      goto getNormalFromElements;
+    }
+  }
+  return;
+
+ getNormalFromElements: ;
+
+/*--------------------------------------------------------------------*
+ * Geometry information cannot be used - generate normals from the
+ * elements
+ *--------------------------------------------------------------------*/
+
+  {
+    if(warnNormFromElem) {
+      Msg::Warning("Some or all boundary normals were determined from mesh "
+                   "elements instead of from the geometry");
+      warnNormFromElem = false;
+    }
+    // Sum the normals from each element connected to the vertex and in the
+    // zone.  Each normal has to be converted independently into an inwards-
+    // pointing normal.
+    //**Weight each normal by the area of the boundary element?
+    SVector3 boNormal(0.);
+    const int nFace = faces.size();
+    for(int iFace = 0; iFace != nFace; ++iFace) {
+      // Normal to the boundary (unknown direction)
+      SVector3 bnt = faces[iFace].parentElement->getFace
+        (faces[iFace].parentFace).normal();
+      // Inwards normal
+      const SVector3 inwards(vertex->point(),
+                             faces[iFace].parentElement->barycenter());
+      if(dot(bnt, inwards) < 0.) bnt.negate();
+      boNormal += bnt;
+    }
+    boNormal.normalize();
+    zoneBoVec.push_back(VertexBoundary(zoneIndex, 0, boNormal,
+                                       const_cast<MVertex*>(vertex),
+                                       vertIndex));
+    patch.add(0);
+  }
+
+}
+
+
+/*******************************************************************************
+ *
+ * Routines from class MZoneBoundary
+ *
+ ******************************************************************************/
+
+
+/*==============================================================================
+ *
+ * Routine: interiorBoundaryVertices
+ *
+ * Purpose
+ * =======
+ *
+ *   Adds a zone to the class.  The boundary vertices from the zone are matched
+ *   against existing vertices in the class and "currently known" connectivity
+ *   for the zone is returned as output.  Records of boundary vertices with
+ *   incomplete connectivity are maintained.
+ *
+ * I/O
+ * ===
+ *
+ *   newZoneIndex       - (I) Index of the new zone (a unique index used to
+ *                            internally label a zone
+ *   mZone              - (I) A description of a zone processed using the
+ *                            MZone class.
+ *   zoneConnMap        - (O) Currently known connectivity for this zone
+ *
+ *============================================================================*/
+
+template <unsigned DIM>
+int MZoneBoundary<DIM>::interiorBoundaryVertices
+(const int newZoneIndex, const MZone<DIM> &mZone, ZoneConnMap &zoneConnMap)
+{
+
+  if(mZone.boVertMap.size() == 0) return 1;
+  zoneConnMap.clear();
+
+  // Loop over all the boundary vertices in 'mZone'
+  for(typename MZone<DIM>::BoVertexMap::const_iterator vMapIt =
+        mZone.boVertMap.begin(); vMapIt != mZone.boVertMap.end(); ++vMapIt) {
+
+//--Find or insert this vertex into the global map
+
+    std::pair<typename GlobalBoVertexMap::iterator, bool> insGblBoVertMap =
+      globalBoVertMap.insert(std::pair<const MVertex*, GlobalVertexData<FaceT> >
+                             (vMapIt->first, GlobalVertexData<FaceT>()));
+    GlobalVertexData<FaceT> &globalVertData = insGblBoVertMap.first->second;
+                                        // Ref. to the global boundary vertex
+                                        // data
+    const ZoneVertexData<typename MZone<DIM>::BoFaceMap::const_iterator>
+      &zoneVertData = vMapIt->second;   // Ref. to boundary vertex data for a
+                                        // zone
+    if(insGblBoVertMap.second) {
+
+//--A new vertex was inserted
+
+      // Copy faces
+      const int nFace = zoneVertData.faces.size();
+      for(int iFace = 0; iFace != nFace; ++iFace) {
+        globalVertData.faces.push_back
+          (typename GlobalVertexData<FaceT>::FaceDataB
+           (newZoneIndex, zoneVertData.faces[iFace]));
+      }
+      // Copy information about the vertex in the zone
+      globalVertData.zoneData.push_back
+        (typename GlobalVertexData<FaceT>::ZoneData
+         (zoneVertData.index, newZoneIndex));
+    }
+    else {
+
+//--This vertex already exists and zone connectivity must be determined
+
+      // Add to the zone connectivity
+      const int nPrevZone = globalVertData.zoneData.size();
+      for(int iPrevZone = 0; iPrevZone != nPrevZone; ++iPrevZone) {
+        ZoneConnectivity &zoneConn =
+          zoneConnMap[ZonePair(globalVertData.zoneData[iPrevZone].zoneIndex,
+                               newZoneIndex)];
+        zoneConn.vertexPairVec.push_back
+          (ZoneConnectivity::VertexPair
+           // const_cast okay, no changes to keys in 'globalBoVertMap'
+           (const_cast<MVertex*>(insGblBoVertMap.first->first),
+            globalVertData.zoneData[iPrevZone].zoneIndex,
+            newZoneIndex,
+            globalVertData.zoneData[iPrevZone].vertexIndex,
+            zoneVertData.index));
+      }
+
+      // Update the list of faces attached to this vertex
+      const int nGFace = globalVertData.faces.size();
+                                        // This is the maximum number of faces
+                                        // searched from 'globalVertData'.  This
+                                        // size may increase as faces are added
+                                        // and the order may also change as
+                                        // faces are deleted.  However, only the
+                                        // faces before nGFace are important.
+      const typename MZone<DIM>::BoFaceMap::const_iterator *zFace =
+        &zoneVertData.faces[0];
+      for(int nZFace = zoneVertData.faces.size(); nZFace--;) {
+        int iGFace = 0;
+        while(iGFace != nGFace) {
+          if((*zFace)->first == globalVertData.faces[iGFace].face) {
+            // Faces match - delete from 'globalVertData'
+            globalVertData.faces.erase(iGFace);
+            break;
+          }
+          ++iGFace;
+        }
+        if(iGFace == nGFace) {
+          // New face - add to 'globalVertData'
+          globalVertData.faces.push_back
+            (typename GlobalVertexData<FaceT>::FaceDataB(newZoneIndex, *zFace));
+        }
+        ++zFace;
+      }
+
+      // If there are no more faces, connectivity for this vertex is complete
+      // and it may be deleted.
+      if(globalVertData.faces.size() == 0)
+        globalBoVertMap.erase(insGblBoVertMap.first);
+    }
+  }  // End loop over boundary vertices
+
+  return 0;
+
+}
+
+
+/*==============================================================================
+ *
+ * Routine: exteriorBoundaryVertices
+ *
+ * Purpose
+ * =======
+ *
+ *   Writes out the remaining unconnected boundary vertices.  These are assumed
+ *   to be at the extent of the domain.
+ *
+ * I/O
+ * ===
+ *
+ *   zoneBoMap          - (O) Boundary vertices for all zones.
+ *
+ * Notes
+ * =====
+ *
+ *   - Should only be called after all zones have been added via
+ *     'interior_boundary_vertices'
+ *
+ ******************************************************************************/
+
+template <unsigned DIM>
+int MZoneBoundary<DIM>::exteriorBoundaryVertices
+(ZoneBoVec &zoneBoVec)
+{
+
+  if(globalBoVertMap.size() == 0) return 1;
+
+  zoneBoVec.clear();
+  zoneBoVec.reserve(3*globalBoVertMap.size()/2);
+
+  BCPatchIndex patch;                   // Provides a BC patch index for each
+                                        // entity
+  bool warnNormFromElem = true;         // A warning that normals were
+                                        // determined from elements.  This
+                                        // warning is only give once (after
+                                        // which the flag is set to false)
+
+  const typename GlobalBoVertexMap::const_iterator vMapEnd =
+    globalBoVertMap.end();
+  for(typename GlobalBoVertexMap::const_iterator vMapIt =
+        globalBoVertMap.begin(); vMapIt != vMapEnd; ++vMapIt) {
+    const int nZone = vMapIt->second.zoneData.size();
+    for(int iZone = 0; iZone != nZone; ++iZone) {
+      const typename GlobalVertexData<FaceT>::ZoneData &zoneData =
+        vMapIt->second.zoneData[iZone]; // Ref. to data stored for this zone
+
+//--Try to find an outwards normal for this vertex
+
+      updateBoVec<DIM, FaceT>(vMapIt->first, zoneData.zoneIndex,
+                              zoneData.vertexIndex, vMapIt->second.faces,
+                              zoneBoVec, patch, warnNormFromElem);
+    }
+  }
+
+  // zoneBoVec has entities as a patch index.  Update with the actual patch
+  // index in 'patch'.
+  patch.generatePatchIndices();
+  const int nBoVert = zoneBoVec.size();
+  for(int iBoVert = 0; iBoVert != nBoVert; ++iBoVert) {
+    zoneBoVec[iBoVert].bcPatchIndex =
+      patch.getIndex(zoneBoVec[iBoVert].bcPatchIndex);
+  }
+
+  return 0;
+
+}
+
+
+/*******************************************************************************
+ *
+ * Explicit instantiations of class MZoneBoundary
+ *
+ ******************************************************************************/
+
+template class MZoneBoundary<2>;
+template class MZoneBoundary<3>;
diff --git a/Geo/MZoneBoundary.h b/Geo/MZoneBoundary.h
new file mode 100644
index 0000000000..9cdf8ba223
--- /dev/null
+++ b/Geo/MZoneBoundary.h
@@ -0,0 +1,266 @@
+// Gmsh - Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
+//
+// See the LICENSE.txt file for license information. Please report all
+// bugs and problems to <gmsh@geuz.org>.
+//
+// MZoneBoundary.h - Copyright (C) 2008 S. Guzik, C. Geuzaine, J.-F. Remacle
+#ifndef _MZONEBOUNDARY_H_
+#define _MZONEBOUNDARY_H_
+
+
+/*******************************************************************************
+ *
+ * - The classes in this file are used to determine connectivity between
+ *   multiple zones and, eventually, the domain boundary vertices
+ *
+ ******************************************************************************/
+
+#include "MZone.h"
+
+
+/*==============================================================================
+ * Required types
+ *============================================================================*/
+
+/*--------------------------------------------------------------------*
+ * Internal zone connectivity
+ *--------------------------------------------------------------------*/
+
+//--Interface between two zones for connectivity
+
+struct ZonePair
+{
+  int zone1;
+  int zone2;
+  ZonePair(const int _zone1, const int _zone2) 
+  {
+    if(_zone1 < zone2) {
+      zone1 = _zone1;
+      zone2 = _zone2;
+    }
+    else {
+      zone1 = _zone2;
+      zone2 = _zone1;
+    }
+  }
+};
+
+inline bool operator==(const ZonePair &zpA, const ZonePair &zpB)
+{
+   return (zpA.zone1 == zpB.zone1 && zpA.zone2 == zpB.zone2);
+}
+
+// Less than for std::map
+struct Less_ZonePair
+  : public std::binary_function<ZonePair, ZonePair, bool> 
+{
+  bool operator()(const ZonePair &zpA, const ZonePair &zpB) const
+  {
+    if(zpA.zone1 < zpB.zone1) return true;
+    if(zpA.zone1 > zpB.zone1) return false;
+    if(zpA.zone2 < zpB.zone2) return true;
+    return false;
+  }
+};
+
+//--Definition of the zone connectivity (a vector of vertex pairs for 2 zones).
+
+struct ZoneConnectivity
+{
+  // Internal structures
+  struct VertexPair                     // Pairs of vertices.  Ordered based on
+                                        // zone indices.
+  {
+    MVertex *vertex;
+    int vertexIndex1;
+    int vertexIndex2;
+    // Constructors
+    VertexPair()
+      : vertex(0), vertexIndex1(0), vertexIndex2(0)
+    { }
+    VertexPair(MVertex *const _vertex,
+               const unsigned zone1, const unsigned zone2,
+               const unsigned _vertexIndex1, const unsigned _vertexIndex2)
+      : vertex(_vertex)
+    {
+      if(zone1 < zone2) {
+        vertexIndex1 = _vertexIndex1;
+        vertexIndex2 = _vertexIndex2;
+      }
+      else {
+        vertexIndex1 = _vertexIndex2;
+        vertexIndex2 = _vertexIndex1;
+      }
+    }
+  };
+  // Data
+  std::vector<VertexPair> vertexPairVec;
+  // Constructor
+  ZoneConnectivity()
+  {
+    vertexPairVec.reserve(32);  // Avoid small reallocations
+  }
+};
+
+//--Output type for zone connectivity
+
+typedef std::map<ZonePair, ZoneConnectivity, Less_ZonePair> ZoneConnMap;
+
+/*--------------------------------------------------------------------*
+ * Boundaries at the domain extent
+ *--------------------------------------------------------------------*/
+
+struct VertexBoundary 
+{
+  int zoneIndex;
+  int bcPatchIndex;
+  SVector3 normal;
+  MVertex *vertex;
+  int vertexIndex;
+  // Constructors
+  VertexBoundary()
+    : vertex(0), vertexIndex(0)
+  { }
+  VertexBoundary(const int _zoneIndex, const int _bcPatchIndex,
+                 const SVector3 &_normal, MVertex *const _vertex,
+                 const int _vertexIndex)
+    : zoneIndex(_zoneIndex), bcPatchIndex(_bcPatchIndex), normal(_normal),
+    vertex(_vertex), vertexIndex(_vertexIndex)
+  { }
+};
+
+typedef std::vector<VertexBoundary> ZoneBoVec;
+
+//--Function object for sorting the ZoneBoVec vector by zone and then BC patch
+//--index
+
+struct ZoneBoVecSort
+{
+  bool operator()(const int i0, const int i1)
+  {
+    if(zoneBoVec[i0].zoneIndex == zoneBoVec[i1].zoneIndex)
+      return zoneBoVec[i0].bcPatchIndex < zoneBoVec[i1].bcPatchIndex;
+    return zoneBoVec[i0].zoneIndex < zoneBoVec[i1].zoneIndex;
+  }
+  ZoneBoVecSort(const ZoneBoVec &_zoneBoVec)
+    : zoneBoVec(_zoneBoVec)
+  { }
+ private:
+  const ZoneBoVec &zoneBoVec;
+};
+
+
+/*******************************************************************************
+ *
+ * class: MZoneBoundary
+ *
+ * Purpose
+ * =======
+ *
+ *   Determines the connectivity between zones (internal boundaries) and
+ *   vertices/elements at the extent of the domain (external boundaries)
+ *
+ *   Template parameters:
+ *     DIM              - dimension of the problem
+ *
+ * Notes
+ * =====
+ *
+ *   - explicitly instantiated in 'MZoneBoundary.cpp'
+ *
+ ******************************************************************************/
+
+template <unsigned DIM>
+class MZoneBoundary
+{
+
+
+/*==============================================================================
+ * Internal types
+ *============================================================================*/
+
+ private:
+
+//--Type of face (MEdge or MFace)
+
+  typedef typename DimTr<DIM>::FaceT FaceT;
+
+//--Data stored for connectivity of vertices
+
+  template<typename FaceT>
+  struct GlobalVertexData
+  {
+    struct FaceDataB
+    {
+      FaceT face;
+      MElement *parentElement;
+      int parentFace;
+      int faceIndex;
+      int zoneIndex;
+      FaceDataB(const int _zoneIndex,
+                const typename MZone<DIM>::BoFaceMap::const_iterator bFMapIt)
+        : face(bFMapIt->first), parentElement(bFMapIt->second.parentElement),
+        parentFace(bFMapIt->second.parentFace),
+        faceIndex(bFMapIt->second.faceIndex), zoneIndex(_zoneIndex)
+      { }
+    };
+    struct ZoneData
+    {
+      int vertexIndex;
+      int zoneIndex;
+      ZoneData(const int _vertexIndex, const int _zoneIndex)
+        : vertexIndex(_vertexIndex), zoneIndex(_zoneIndex)
+      { }
+    };
+    CCon::FaceVector<FaceDataB> faces;
+    CCon::FaceVector<ZoneData> zoneData;
+                                        // A 'FaceVector' is not strictly
+                                        // optimized for the vertices but should
+                                        // still work quite well.
+    // Constructor
+    GlobalVertexData() { }
+  };
+
+  typedef std::map<const MVertex*, GlobalVertexData<FaceT>,
+    std::less<const MVertex*> > GlobalBoVertexMap;
+
+
+/*==============================================================================
+ * Member functions
+ *============================================================================*/
+
+ public:
+
+//--Reset the database
+
+  void clear()
+  {
+    globalBoVertMap.clear();
+  }
+
+//--Add a zone to the global map of boundary vertices and return connectivity
+//--between zones.
+
+  int interiorBoundaryVertices(const int newZoneIndex, const MZone<DIM> &mZone,
+                               ZoneConnMap &zoneConnMap);
+  
+//--Return exterior boundary vertices (unconnected vertices at the extent of the
+//--domain
+
+  int exteriorBoundaryVertices(ZoneBoVec &zoneBoVec);
+
+
+/*==============================================================================
+ * Member data
+ *============================================================================*/
+
+private:
+
+//--Data members
+
+  GlobalBoVertexMap globalBoVertMap;    // Map of unconnected boundary vertices
+                                        // for the entire domain
+
+};
+
+#endif
diff --git a/Geo/Makefile b/Geo/Makefile
index 642df5f802..8e9ef49de6 100644
--- a/Geo/Makefile
+++ b/Geo/Makefile
@@ -38,7 +38,8 @@ SRC = GEntity.cpp\
       GaussQuadratureHex.cpp\
       GaussLegendreSimplex.cpp\
       MFace.cpp\
-      MElement.cpp
+      MElement.cpp\
+      MZone.cpp MZoneBoundary.cpp
 
 OBJ = ${SRC:.cpp=${OBJEXT}}
 
@@ -68,7 +69,7 @@ depend:
 # DO NOT DELETE THIS LINE
 GEntity.o: GEntity.cpp GEntity.h Range.h SPoint3.h SBoundingBox3d.h \
   ../Common/VertexArray.h ../Geo/SVector3.h ../Geo/SPoint3.h \
-  ../Common/Context.h
+  ../Common/Context.h ../Geo/CGNSOptions.h
 GVertex.o: GVertex.cpp GVertex.h GEntity.h Range.h SPoint3.h \
   SBoundingBox3d.h GPoint.h SPoint2.h GFace.h GEdgeLoop.h GEdge.h \
   SVector3.h Pair.h MVertex.h ../Common/Message.h
@@ -85,7 +86,8 @@ GFace.o: GFace.cpp GModel.h GVertex.h GEntity.h Range.h SPoint3.h \
   GEdgeLoop.h Pair.h GRegion.h MElement.h ../Common/GmshDefines.h \
   MVertex.h MEdge.h MFace.h ../Common/Message.h ../Numeric/Numeric.h \
   ../Numeric/NumericEmbedded.h ../Numeric/GaussLegendre1D.h \
-  ../Common/VertexArray.h ../Geo/SVector3.h ../Common/Context.h
+  ../Common/VertexArray.h ../Geo/SVector3.h ../Common/Context.h \
+  ../Geo/CGNSOptions.h
 GRegion.o: GRegion.cpp GModel.h GVertex.h GEntity.h Range.h SPoint3.h \
   SBoundingBox3d.h GPoint.h SPoint2.h GEdge.h SVector3.h GFace.h \
   GEdgeLoop.h Pair.h GRegion.h MElement.h ../Common/GmshDefines.h \
@@ -103,14 +105,16 @@ gmshEdge.o: gmshEdge.cpp GModel.h GVertex.h GEntity.h Range.h SPoint3.h \
   gmshSurface.h ../Numeric/Numeric.h ../Numeric/NumericEmbedded.h \
   ../Common/ListUtils.h ../Common/TreeUtils.h ../Common/avl.h \
   ../Common/ListUtils.h ExtrudeParams.h ../Common/SmoothData.h \
-  GeoInterpolation.h ../Common/Message.h ../Common/Context.h
+  GeoInterpolation.h ../Common/Message.h ../Common/Context.h \
+  ../Geo/CGNSOptions.h
 gmshFace.o: gmshFace.cpp GModel.h GVertex.h GEntity.h Range.h SPoint3.h \
   SBoundingBox3d.h GPoint.h SPoint2.h GEdge.h SVector3.h GFace.h \
   GEdgeLoop.h Pair.h GRegion.h gmshFace.h Geo.h ../Common/GmshDefines.h \
   gmshSurface.h ../Numeric/Numeric.h ../Numeric/NumericEmbedded.h \
   ../Common/ListUtils.h ../Common/TreeUtils.h ../Common/avl.h \
   ../Common/ListUtils.h ExtrudeParams.h ../Common/SmoothData.h \
-  GeoInterpolation.h ../Common/Message.h ../Common/Context.h
+  GeoInterpolation.h ../Common/Message.h ../Common/Context.h \
+  ../Geo/CGNSOptions.h
 gmshRegion.o: gmshRegion.cpp GModel.h GVertex.h GEntity.h Range.h \
   SPoint3.h SBoundingBox3d.h GPoint.h SPoint2.h GEdge.h SVector3.h \
   GFace.h GEdgeLoop.h Pair.h GRegion.h gmshRegion.h Geo.h \
@@ -120,7 +124,8 @@ gmshRegion.o: gmshRegion.cpp GModel.h GVertex.h GEntity.h Range.h \
   ExtrudeParams.h ../Common/SmoothData.h ../Common/Message.h
 gmshSurface.o: gmshSurface.cpp gmshSurface.h Pair.h Range.h SPoint2.h \
   SPoint3.h SVector3.h SBoundingBox3d.h ../Numeric/Numeric.h \
-  ../Numeric/NumericEmbedded.h ../Common/Message.h
+  ../Numeric/NumericEmbedded.h ../Common/Message.h \
+  ../contrib/MathEval/matheval.h
 OCCVertex.o: OCCVertex.cpp GModel.h GVertex.h GEntity.h Range.h SPoint3.h \
   SBoundingBox3d.h GPoint.h SPoint2.h GEdge.h SVector3.h GFace.h \
   GEdgeLoop.h Pair.h GRegion.h OCCVertex.h OCCIncludes.h MVertex.h \
@@ -128,13 +133,14 @@ OCCVertex.o: OCCVertex.cpp GModel.h GVertex.h GEntity.h Range.h SPoint3.h \
 OCCEdge.o: OCCEdge.cpp GModel.h GVertex.h GEntity.h Range.h SPoint3.h \
   SBoundingBox3d.h GPoint.h SPoint2.h GEdge.h SVector3.h GFace.h \
   GEdgeLoop.h Pair.h GRegion.h ../Common/Message.h OCCEdge.h OCCVertex.h \
-  OCCIncludes.h MVertex.h OCCFace.h ../Common/Context.h
+  OCCIncludes.h MVertex.h OCCFace.h ../Common/Context.h \
+  ../Geo/CGNSOptions.h
 OCCFace.o: OCCFace.cpp GModel.h GVertex.h GEntity.h Range.h SPoint3.h \
   SBoundingBox3d.h GPoint.h SPoint2.h GEdge.h SVector3.h GFace.h \
   GEdgeLoop.h Pair.h GRegion.h OCCVertex.h OCCIncludes.h MVertex.h \
   OCCEdge.h OCCFace.h ../Common/Message.h ../Numeric/Numeric.h \
   ../Numeric/NumericEmbedded.h ../Common/VertexArray.h ../Geo/SVector3.h \
-  ../Common/Context.h
+  ../Common/Context.h ../Geo/CGNSOptions.h
 OCCRegion.o: OCCRegion.cpp GModel.h GVertex.h GEntity.h Range.h SPoint3.h \
   SBoundingBox3d.h GPoint.h SPoint2.h GEdge.h SVector3.h GFace.h \
   GEdgeLoop.h Pair.h GRegion.h OCCVertex.h OCCIncludes.h MVertex.h \
@@ -163,7 +169,7 @@ discreteRegion.o: discreteRegion.cpp discreteRegion.h GModel.h GVertex.h \
 fourierEdge.o: fourierEdge.cpp fourierEdge.h GEdge.h GEntity.h Range.h \
   SPoint3.h SBoundingBox3d.h GVertex.h GPoint.h SPoint2.h SVector3.h \
   GModel.h GFace.h GEdgeLoop.h Pair.h GRegion.h fourierVertex.h MVertex.h \
-  ../Common/Context.h
+  ../Common/Context.h ../Geo/CGNSOptions.h
 fourierFace.o: fourierFace.cpp fourierVertex.h GModel.h GVertex.h \
   GEntity.h Range.h SPoint3.h SBoundingBox3d.h GPoint.h SPoint2.h GEdge.h \
   SVector3.h GFace.h GEdgeLoop.h Pair.h GRegion.h MVertex.h fourierFace.h \
@@ -172,7 +178,7 @@ fourierProjectionFace.o: fourierProjectionFace.cpp \
   fourierProjectionFace.h GModel.h GVertex.h GEntity.h Range.h SPoint3.h \
   SBoundingBox3d.h GPoint.h SPoint2.h GEdge.h SVector3.h GFace.h \
   GEdgeLoop.h Pair.h GRegion.h ../Common/VertexArray.h ../Geo/SVector3.h \
-  ../Common/Context.h
+  ../Common/Context.h ../Geo/CGNSOptions.h
 GModel.o: GModel.cpp GModel.h GVertex.h GEntity.h Range.h SPoint3.h \
   SBoundingBox3d.h GPoint.h SPoint2.h GEdge.h SVector3.h GFace.h \
   GEdgeLoop.h Pair.h GRegion.h MElement.h ../Common/GmshDefines.h \
@@ -183,7 +189,8 @@ GModel.o: GModel.cpp GModel.h GVertex.h GEntity.h Range.h SPoint3.h \
   ../Geo/gmshSurface.h ../Common/ListUtils.h ../Common/TreeUtils.h \
   ../Common/avl.h ../Common/ListUtils.h ../Geo/SPoint2.h \
   ../Geo/ExtrudeParams.h ../Common/SmoothData.h ../Post/PView.h \
-  ../Geo/SPoint3.h ../Mesh/Generator.h ../Common/Context.h
+  ../Geo/SPoint3.h ../Mesh/Generator.h ../Common/Context.h \
+  ../Geo/CGNSOptions.h
 GModelIO_Geo.o: GModelIO_Geo.cpp GModel.h GVertex.h GEntity.h Range.h \
   SPoint3.h SBoundingBox3d.h GPoint.h SPoint2.h GEdge.h SVector3.h \
   GFace.h GEdgeLoop.h Pair.h GRegion.h Geo.h ../Common/GmshDefines.h \
@@ -202,9 +209,9 @@ GModelIO_Mesh.o: GModelIO_Mesh.cpp GModel.h GVertex.h GEntity.h Range.h \
 GModelIO_OCC.o: GModelIO_OCC.cpp GModelIO_OCC.h GModel.h GVertex.h \
   GEntity.h Range.h SPoint3.h SBoundingBox3d.h GPoint.h SPoint2.h GEdge.h \
   SVector3.h GFace.h GEdgeLoop.h Pair.h GRegion.h OCCIncludes.h \
-  ../Common/Message.h ../Common/Context.h OCCVertex.h MVertex.h OCCEdge.h \
-  OCCFace.h OCCRegion.h MElement.h ../Common/GmshDefines.h MEdge.h \
-  MFace.h ../Common/OpenFile.h
+  ../Common/Message.h ../Common/Context.h ../Geo/CGNSOptions.h \
+  OCCVertex.h MVertex.h OCCEdge.h OCCFace.h OCCRegion.h MElement.h \
+  ../Common/GmshDefines.h MEdge.h MFace.h ../Common/OpenFile.h
 GModelIO_Fourier.o: GModelIO_Fourier.cpp GModel.h GVertex.h GEntity.h \
   Range.h SPoint3.h SBoundingBox3d.h GPoint.h SPoint2.h GEdge.h \
   SVector3.h GFace.h GEdgeLoop.h Pair.h GRegion.h ../Common/Message.h \
@@ -212,9 +219,10 @@ GModelIO_Fourier.o: GModelIO_Fourier.cpp GModel.h GVertex.h GEntity.h \
   GModelIO_Fourier.h
 GModelIO_CGNS.o: GModelIO_CGNS.cpp GModel.h GVertex.h GEntity.h Range.h \
   SPoint3.h SBoundingBox3d.h GPoint.h SPoint2.h GEdge.h SVector3.h \
-  GFace.h GEdgeLoop.h Pair.h GRegion.h ../Common/Message.h MElement.h \
-  ../Common/GmshDefines.h MVertex.h MEdge.h MFace.h MNeighbour.h \
-  MEdgeHash.h ../Common/Hash.h MFaceHash.h
+  GFace.h GEdgeLoop.h Pair.h GRegion.h ../Common/Message.h MZone.h \
+  MElement.h ../Common/GmshDefines.h MVertex.h MEdge.h MFace.h \
+  MEdgeHash.h ../Common/Hash.h MFaceHash.h CustomContainer.h \
+  MZoneBoundary.h CGNSOptions.h
 GModelIO_MED.o: GModelIO_MED.cpp GModel.h GVertex.h GEntity.h Range.h \
   SPoint3.h SBoundingBox3d.h GPoint.h SPoint2.h GEdge.h SVector3.h \
   GFace.h GEdgeLoop.h Pair.h GRegion.h ../Common/Message.h
@@ -232,7 +240,7 @@ Geo.o: Geo.cpp ../Common/Message.h ../Numeric/Numeric.h \
   ExtrudeParams.h ../Common/SmoothData.h GModel.h GVertex.h GEntity.h \
   GPoint.h GEdge.h GFace.h GEdgeLoop.h GRegion.h GeoInterpolation.h \
   ../Mesh/Field.h ../Geo/Geo.h ../Post/PView.h ../Geo/SPoint3.h \
-  ../Common/Context.h
+  ../Common/Context.h ../Geo/CGNSOptions.h
 GeoStringInterface.o: GeoStringInterface.cpp ../Common/Message.h \
   ../Numeric/Numeric.h ../Numeric/NumericEmbedded.h \
   ../Common/MallocUtils.h ../Common/StringUtils.h Geo.h \
@@ -240,8 +248,9 @@ GeoStringInterface.o: GeoStringInterface.cpp ../Common/Message.h \
   SPoint3.h SVector3.h SBoundingBox3d.h ../Common/ListUtils.h \
   ../Common/TreeUtils.h ../Common/avl.h ../Common/ListUtils.h \
   ExtrudeParams.h ../Common/SmoothData.h GeoStringInterface.h \
-  ../Common/OpenFile.h ../Common/Context.h GModel.h GVertex.h GEntity.h \
-  GPoint.h GEdge.h GFace.h GEdgeLoop.h GRegion.h ../Parser/Parser.h
+  ../Common/OpenFile.h ../Common/Context.h ../Geo/CGNSOptions.h GModel.h \
+  GVertex.h GEntity.h GPoint.h GEdge.h GFace.h GEdgeLoop.h GRegion.h \
+  ../Parser/Parser.h
 GeoInterpolation.o: GeoInterpolation.cpp ../Common/Message.h Geo.h \
   ../Common/GmshDefines.h gmshSurface.h Pair.h Range.h SPoint2.h \
   SPoint3.h SVector3.h SBoundingBox3d.h ../Numeric/Numeric.h \
@@ -252,7 +261,7 @@ GeoInterpolation.o: GeoInterpolation.cpp ../Common/Message.h Geo.h \
 findLinks.o: findLinks.cpp ../Common/Message.h GModel.h GVertex.h \
   GEntity.h Range.h SPoint3.h SBoundingBox3d.h GPoint.h SPoint2.h GEdge.h \
   SVector3.h GFace.h GEdgeLoop.h Pair.h GRegion.h ../Common/TreeUtils.h \
-  ../Common/avl.h ../Common/ListUtils.h
+  ../Common/avl.h ../Common/ListUtils.h ../Common/ListUtils.h
 MVertex.o: MVertex.cpp MVertex.h SPoint3.h GEdge.h GEntity.h Range.h \
   SBoundingBox3d.h GVertex.h GPoint.h SPoint2.h SVector3.h GFace.h \
   GEdgeLoop.h Pair.h ../Common/Message.h
@@ -272,11 +281,22 @@ GaussLegendreSimplex.o: GaussLegendreSimplex.cpp MElement.h \
   ../Common/GmshDefines.h MVertex.h SPoint3.h MEdge.h SVector3.h MFace.h \
   GaussLegendreSimplex.h ../Numeric/GaussLegendre1D.h
 MFace.o: MFace.cpp MFace.h MVertex.h SPoint3.h SVector3.h \
-  ../Numeric/Numeric.h ../Numeric/NumericEmbedded.h ../Common/Context.h
+  ../Numeric/Numeric.h ../Numeric/NumericEmbedded.h ../Common/Context.h \
+  ../Geo/CGNSOptions.h
 MElement.o: MElement.cpp MElement.h ../Common/GmshDefines.h MVertex.h \
   SPoint3.h MEdge.h SVector3.h MFace.h GEntity.h Range.h SBoundingBox3d.h \
   GFace.h GPoint.h GEdgeLoop.h GEdge.h GVertex.h SPoint2.h Pair.h \
   ../Numeric/Numeric.h ../Numeric/NumericEmbedded.h \
   ../Numeric/FunctionSpace.h ../Common/GmshMatrix.h \
   ../Numeric/GaussLegendre1D.h ../Common/Message.h ../Common/Context.h \
-  ../Mesh/qualityMeasures.h
+  ../Geo/CGNSOptions.h ../Mesh/qualityMeasures.h
+MZone.o: MZone.cpp MZone.h MElement.h ../Common/GmshDefines.h MVertex.h \
+  SPoint3.h MEdge.h SVector3.h MFace.h GFace.h GEntity.h Range.h \
+  SBoundingBox3d.h GPoint.h GEdgeLoop.h GEdge.h GVertex.h SPoint2.h \
+  Pair.h GRegion.h MEdgeHash.h ../Common/Hash.h MFaceHash.h \
+  CustomContainer.h
+MZoneBoundary.o: MZoneBoundary.cpp MZoneBoundary.h MZone.h MElement.h \
+  ../Common/GmshDefines.h MVertex.h SPoint3.h MEdge.h SVector3.h MFace.h \
+  GFace.h GEntity.h Range.h SBoundingBox3d.h GPoint.h GEdgeLoop.h GEdge.h \
+  GVertex.h SPoint2.h Pair.h GRegion.h MEdgeHash.h ../Common/Hash.h \
+  MFaceHash.h CustomContainer.h ../Common/Message.h
diff --git a/Geo/SVector3.h b/Geo/SVector3.h
index 5de09ab41e..f3f0d90637 100644
--- a/Geo/SVector3.h
+++ b/Geo/SVector3.h
@@ -30,6 +30,7 @@ class SVector3 {
     double n = norm(); if(n){ P[0] /= n; P[1] /= n; P[2] /= n; }
     return n;
   }
+  void negate() { P[0] = -P[0]; P[1] = -P[1]; P[2] = -P[2]; }  
   // why both [] and (), why not
   double &operator[](int i){ return P[i]; }
   double operator[](int i) const { return P[i]; }
@@ -37,20 +38,22 @@ class SVector3 {
   double operator()(int i) const { return P[i]; }
   SVector3 & operator += (const SVector3 &a)
   {
-    for(int i = 0; i < 3; i++)
-      P[i] += a[i];
+    P[0] += a[0];  P[1] += a[1];  P[2] += a[2];
     return *this;
   }
   SVector3 & operator -= (const SVector3 &a)
   { 
-    for(int i = 0; i < 3; i++)
-      P[i] -= a[i];
+    P[0] -= a[0];  P[1] -= a[1];  P[2] -= a[2];
     return *this;
   }
   SVector3 & operator *= (const SVector3 &a)
-  { 
-    for(int i = 0; i < 3; i++)
-      P[i] *= a[i];
+  {
+    P[0] *= a[0]; P[1] *= a[1]; P[2] *= a[2];
+    return *this;
+  }
+  SVector3 & operator *= (const double v)
+  {
+    P[0] *= v; P[1] *= v; P[2] *= v;
     return *this;
   }
   SVector3 & operator = (double v)
-- 
GitLab