diff --git a/Geo/GModelIO_CGNS.cpp b/Geo/GModelIO_CGNS.cpp
index 0862f2a7126e608f2be9ded7478c0b6ca6c65fda..4ed7f00a551ab0e4110c738d1497149885afccbe 100644
--- a/Geo/GModelIO_CGNS.cpp
+++ b/Geo/GModelIO_CGNS.cpp
@@ -399,13 +399,13 @@ int GModel::writeCGNS(const std::string &name, int zoneDefinition,
   if(meshDim == 2) vectorDim = options.vectorDim;
 
   // open the file
-  int cgIndexFile;
+  int cgIndexFile=0;
   if(cg_open(name.c_str(), MODE_WRITE, &cgIndexFile)) return cgnsErr();
 
   // write the base node
-  int cgIndexBase;
-  if(cg_base_write(cgIndexFile, options.baseName.c_str(), meshDim, meshDim,
-                   &cgIndexBase))
+  int cgIndexBase=0;
+  if(cg_base_write(cgIndexFile, options.baseName.c_str(), meshDim, meshDim, 
+                   &cgIndexBase)) 
     return cgnsErr();
 
   // write information about who generated the mesh
@@ -415,17 +415,21 @@ int GModel::writeCGNS(const std::string &name, int zoneDefinition,
   switch(meshDim) {
   case 2:
     MZone<2>::preInit();
-    MZoneBoundary<2>::preInit(); // Add?:
-    write_CGNS_zones<2>(*this, zoneDefinition, ...);
+    MZoneBoundary<2>::preInit();
+    write_CGNS_zones<2>(*this, zoneDefinition, numZone, options, 
+                        scalingFactor, vectorDim, groups[face], 
+                        cgIndexFile, cgIndexBase);
     MZone<2>::postDestroy();
-    MZoneBoundary<2>::postDestroy(); // Add?:
+    MZoneBoundary<2>::postDestroy();
     break;
   case 3:
     MZone<3>::preInit();
-    MZoneBoundary<3>::preInit(); // Add?:
-    write_CGNS_zones<3>(*this, zoneDefinition, ...);
+    MZoneBoundary<3>::preInit();
+    write_CGNS_zones<3>(*this, zoneDefinition, numZone, options, 
+                        scalingFactor, vectorDim, groups[region], 
+                        cgIndexFile, cgIndexBase);
     MZone<3>::postDestroy();
-    MZoneBoundary<3>::postDestroy(); // Add?:
+    MZoneBoundary<3>::postDestroy();
     break;
   }
 
@@ -645,6 +649,8 @@ int get_zone_definition(GModel &model, const int zoneDefinition,
 
   int status = 0;
   const char *_zoneName = "Partition";
+  // NBN: pass name to expand_name(): Zone_&N& ==> Zone_Fluid
+  std::string zn;
 
 //--Get indices for the zonex
 
@@ -665,8 +671,9 @@ int get_zone_definition(GModel &model, const int zoneDefinition,
         break;
       case 2:  // Zone defined by physical
         partition = -1;
-        _zoneName = model.getPhysicalName(meshDim, globalPhysicalIt->first)
-           .c_str();
+      //_zoneName = model.getPhysicalName(meshDim, globalPhysicalIt->first).c_str();
+        zn = model.getPhysicalName(meshDim, globalPhysicalIt->first);   // NBN:
+        _zoneName = zn.c_str();   // NBN:
         physicalItBegin = globalPhysicalIt++;
         physicalItEnd = globalPhysicalIt;
         break;
@@ -809,8 +816,8 @@ int write_CGNS_zones(GModel &model, const int zoneDefinition, const int numZone,
 //--Write the zone
 
           // Write the zone node
-          int cgIndexZone;
-          int cgZoneSize[3];
+          int cgIndexZone=0;
+          int cgZoneSize[3]={0};
           cgZoneSize[0] = writeZone->zoneVertVec.size();  // Number of vertices
 #ifdef CGNS_TEST1
           // Count all the sub-elements in a Triangle 10
@@ -832,40 +839,49 @@ int write_CGNS_zones(GModel &model, const int zoneDefinition, const int numZone,
           if(cg_zone_write(cgIndexFile, cgIndexBase,
                            writeTask->zoneName.c_str(), cgZoneSize,
                            Unstructured, &cgIndexZone))
+          {
             return cgnsErr();
+          }
           // 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));
+          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;
+          int cgIndexGrid=0;
           if(cg_grid_write(cgIndexFile, cgIndexBase, cgIndexZone,
                            "GridCoordinates", &cgIndexGrid))
             return cgnsErr();
 
           // Write the grid coordinates
-          int cgIndexCoord;
+          int cgIndexCoord=0;
           dBuffer.resize(cgZoneSize[0]);
-          // x
-          for(int i = 0; i != cgZoneSize[0]; ++i) 
-              dBuffer[i] = writeZone->zoneVertVec[i]->x()*scalingFactor;
+
+          // x-coordinates for this zone
+          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)
+
+          // y-coordinates for this zone
+          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
+
+          // z-coordinates for this zone
           if(vectorDim == 3) {
-            for(int i = 0; i != cgZoneSize[0]; ++i)
+            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();
@@ -951,7 +967,9 @@ int write_CGNS_zones(GModel &model, const int zoneDefinition, const int numZone,
                   writeZone->zoneElemConn[typeMSHm1].numBoElem + iElemSection,
                   &writeZone->zoneElemConn[typeMSHm1].connectivity[0],
                   &cgIndexSection))
+              {
                 return cgnsErr();
+              }
               ++iElemSection;
             }
           }
@@ -963,7 +981,8 @@ int write_CGNS_zones(GModel &model, const int zoneDefinition, const int numZone,
           // Write the connectivity for each zone pair
           const ZoneConnMap::const_iterator gCEnd = zoneConnMap.end();
           for(ZoneConnMap::const_iterator gCIt = zoneConnMap.begin();
-              gCIt != gCEnd; ++gCIt) {
+              gCIt != gCEnd; ++gCIt) 
+          {
             const int nVert = gCIt->second.vertexPairVec.size();
             iBuffer1.resize(nVert);
             iBuffer2.resize(nVert);
@@ -987,7 +1006,9 @@ int write_CGNS_zones(GModel &model, const int zoneDefinition, const int numZone,
                 &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,
@@ -995,7 +1016,9 @@ int write_CGNS_zones(GModel &model, const int zoneDefinition, const int numZone,
                 &iBuffer2[0], zoneInfo[gCIt->first.zone1].name.c_str(),
                 Unstructured, PointListDonor, Integer, nVert, &iBuffer1[0],
                 &cgIndexInterface))
+            {
               return cgnsErr();
+            }
           }
 
 //--Pop from the queue
@@ -1021,7 +1044,8 @@ int write_CGNS_zones(GModel &model, const int zoneDefinition, const int numZone,
           if(get_zone_definition
              (model, zoneDefinition, numZone, options, DIM, group,
               globalZoneIndex, globalPhysicalIt, zoneTask.zoneIndex, partition,
-              physicalItBegin, physicalItEnd, zoneTask.zoneName)) {
+              physicalItBegin, physicalItEnd, zoneTask.zoneName)) 
+          {
             omp_set_lock(&threadWLock);
             --threadsWorking;
             omp_unset_lock(&threadWLock);
@@ -1030,7 +1054,8 @@ int write_CGNS_zones(GModel &model, const int zoneDefinition, const int numZone,
             zoneTask.zone.clear();
             // Loop through all physical in the zone definition
             for(PhysGroupMap::const_iterator physicalIt = physicalItBegin;
-                physicalIt != physicalItEnd; ++physicalIt) {
+                physicalIt != physicalItEnd; ++physicalIt) 
+            {
               // Add elements from all entities in the physical with defined
               // partition number
               if(partition == -1) {
@@ -1061,67 +1086,80 @@ int write_CGNS_zones(GModel &model, const int zoneDefinition, const int numZone,
  * Write the remaining unconnected vertices as boundary conditions
  *--------------------------------------------------------------------*/
 
-      if(options.writeBC) {
+      if(options.writeBC) 
+      {
         ZoneBoVec zoneBoVec;            // from 'MZoneBoundary.h'
         if(mZoneBoundary.exteriorBoundaryVertices
-           (options.normalSource, zoneBoVec) == 0) {
-
+           (options.normalSource, zoneBoVec) == 0) 
+        {
           // 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-1],
-                                         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;
-            std::string patchName = options.patchName;
-            expand_name(patchName, patchIndex, "Patch");
-            if(patchName.length() == 0) {
-              patchName = "Patch_";
-              patchName += CGNSNameStr(patchIndex+1).c_str();
-            }
-            if(cg_boco_write(cgIndexFile, cgIndexBase,
-                             zoneInfo[zoneIndex].cgIndex, patchName.c_str(),
-                             BCTypeNull, PointList, numBCVert, &iBuffer1[0],
-                             &cgIndexBoco))
-              return cgnsErr();
-
-            if(options.normalSource) {
-              int normalIndex;
-              if(cg_boco_normal_write(cgIndexFile, cgIndexBase,
-                                      zoneInfo[zoneIndex].cgIndex, cgIndexBoco,
-                                      &normalIndex, 1, RealDouble, &dBuffer[0]))
+          if (numBoVert>=1) 
+          {
+            Msg::Info("writing BoVerts...");
+
+            std::vector<int> iZBV(numBoVert);
+            for(int i = 0; i != numBoVert; ++i) iZBV[i] = i;
+            std::sort<int*, ZoneBoVecSort>(&iZBV[0], &iZBV[numBoVert-1],
+                                           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;
+              std::string patchName = options.patchName;
+              expand_name(patchName, patchIndex, "Patch");
+              if(patchName.length() == 0) {
+                patchName = "Patch_";
+                patchName += CGNSNameStr(patchIndex+1).c_str();
+              }
+              if(cg_boco_write(cgIndexFile, cgIndexBase,
+                               zoneInfo[zoneIndex].cgIndex, patchName.c_str(),
+                               BCTypeNull, PointList, numBCVert, &iBuffer1[0],
+                               &cgIndexBoco))
+              {
                 return cgnsErr();
+              }
+
+              if(options.normalSource) {
+                int normalIndex;
+                if(cg_boco_normal_write(cgIndexFile, cgIndexBase,
+                                        zoneInfo[zoneIndex].cgIndex, cgIndexBoco,
+                                        &normalIndex, 1, RealDouble, &dBuffer[0]))
+                {
+                  return cgnsErr();
+                }
+              }
             }
           }
         }
       }
 
-//       std::cout << "Leaving master thread\n";  // DBG
-    }  // End master thread instructions
+      // NBN: MZoneBoundary defines no destructor, so force 
+      // the deallocation of the new pointers with clear()
+      mZoneBoundary.clear();
+      Msg::Info("Leaving master thread");  // DBG
 
+    }  // End master thread instructions
   }  // End omp parallel section
 
 //--Destroy omp locks
@@ -1129,6 +1167,7 @@ int write_CGNS_zones(GModel &model, const int zoneDefinition, const int numZone,
   omp_destroy_lock(&threadWLock);
   omp_destroy_lock(&queueLock);
 
+  return 0;
 }
 
 #else
diff --git a/Geo/MZoneBoundary.cpp b/Geo/MZoneBoundary.cpp
index b5cee4541c28ba66ad02cf22a78a211044500703..0e66c86c95daade8c814bb0134a7c57305698b78 100644
--- a/Geo/MZoneBoundary.cpp
+++ b/Geo/MZoneBoundary.cpp
@@ -9,6 +9,8 @@
 
 #if defined(HAVE_LIBCGNS)
 
+#define DEBUG_FaceT 0       // NBN: toggle reporting of face insertions
+
 #include <iostream> // DBG
 #include <limits> // ?
 
@@ -217,7 +219,7 @@ int edge_normal
  &faces, SVector3 &boNormal, const int onlyFace = -1)
 {
 
-  double par;
+  double par=0.0;
   // Note: const_cast used to match MVertex.cpp interface
   if(!reparamMeshVertexOnEdge(vertex, gEdge, par)) return 1;
 
@@ -580,7 +582,7 @@ void updateBoVec<3, MFace>
             //   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) {
+            if(!matchedFace && (3 == nVOnF)) {
               const MVertex *vertex2;
               const MVertex *vertex3;
               switch(vertexOnF) {
@@ -811,25 +813,54 @@ int MZoneBoundary<DIM>::interiorBoundaryVertices
 
 //--A new vertex was inserted
 
-//       if(debug) {
-//         std::cout << "This vertex is new and has bo faces:\n";
-//       }
       // Copy faces
       const int nFace = zoneVertData.faces.size();
+      int iCount=0;
       for(int iFace = 0; iFace != nFace; ++iFace) {
-//         if(debug) {
-//           std::cout << "  ("
-//                     << zoneVertData.faces[iFace]->first.getVertex(0)->x() << ","
-//                     << zoneVertData.faces[iFace]->first.getVertex(0)->y()
-//                     << "), ("
-//                     << zoneVertData.faces[iFace]->first.getVertex(1)->x() << ","
-//                     << zoneVertData.faces[iFace]->first.getVertex(1)->y()
-//                     << ")\n";
-//         }
+
+#if (0)
+        // NBN: changing the member object for a pointer means 
+        // we need to do the following in two steps (see below)
+
         globalVertData.faces.push_back
           (typename GlobalVertexData<FaceT>::FaceDataB
            (newZoneIndex, zoneVertData.faces[iFace]));
+
+#else
+        // Using FaceAllocator<T> with face2Pool, the std constructors 
+        // are not called, so FaceDataB std members are incomplete.
+
+        // By replacing the FaceT member with FaceT* we fix the issue of 
+        // auto-deleting "un-constructed" containers.  Here, the simple 
+        // data type (pointer) member is allocated in the ctor, then we
+        // create and store its FaceT object once the FaceDataB object 
+        // is safely in the container.
+        //
+        // Note: we must now make sure to delete these pointers.
+        // See adjusted version of MZoneBoundary::clear();
+
+        // Step 1: append new FaceDataB<> object
+        typename GlobalVertexData<FaceT>::FaceDataB& tFDB = 
+          globalVertData.faces.push_back
+              (typename GlobalVertexData<FaceT>::FaceDataB
+              (newZoneIndex, zoneVertData.faces[iFace]) );
+        
+        // Step 2: construct its internal face object
+        tFDB.face = new FaceT(zoneVertData.faces[iFace]->first);
+
+      #if (DEBUG_FaceT)
+        if (! tFDB.face) {
+          Msg::Info("MZoneBoundary<DIM>::interiorBoundaryVertices, failed to alloc face object");
+        } else {
+          int nv = tFDB.face->getNumVertices();
+          Msg::Info("interiorBoundaryVertices: allocated FaceT object %5d with %d verts", ++iCount, nv);
+        }
+      #endif
+
+
+#endif
       }
+
       // Copy information about the vertex in the zone
       globalVertData.zoneData.push_back
         (typename GlobalVertexData<FaceT>::ZoneData
@@ -871,8 +902,11 @@ int MZoneBoundary<DIM>::interiorBoundaryVertices
         &zoneVertData.faces[0];
       for(int nZFace = zoneVertData.faces.size(); nZFace--;) {
         bool foundMatch = false;
-        for(int iGFace = 0; iGFace != nGFace; ++iGFace) {
-          if((*zFace)->first == globalVertData.faces[iGFace].face) {
+        for(int iGFace = 0; iGFace != nGFace; ++iGFace) 
+        {
+          // NBN: face is now a pointer, so need to de-reference
+        //if((*zFace)->first ==  globalVertData.faces[iGFace].face ) 
+          if((*zFace)->first == *(globalVertData.faces[iGFace].face)) {
             foundMatch = true;
             // Faces match - delete from 'globalVertData'.
             globalVertData.faces.erase(iGFace);
@@ -892,8 +926,9 @@ int MZoneBoundary<DIM>::interiorBoundaryVertices
 
       // If there are no more faces, connectivity for this vertex is complete
       // and it may be deleted.
-      if(globalVertData.faces.size() == 0)
+      if(globalVertData.faces.size() == 0) {
         globalBoVertMap.erase(insGblBoVertMap.first);
+      }
       else {
         // Update the list of zones attached to this vertex
         globalVertData.zoneData.push_back
@@ -1012,13 +1047,25 @@ int MZoneBoundary<DIM>::exteriorBoundaryVertices
 template<>
 template<>
 MZoneBoundary<2>::GlobalVertexData<MEdge>::FaceDataB::FaceDataB()
-  : face(0, 0) 
+  : 
+  //face(0, 0),         // NBN: replaced this MEdge object
+    face(NULL),         // NBN: with a pointer to MEdge
+    parentElement(0),   // NBN: also, init members
+    parentFace(0),
+    faceIndex(0),
+    zoneIndex(0)
 { }
 
 template<>
 template<>
 MZoneBoundary<3>::GlobalVertexData<MFace>::FaceDataB::FaceDataB()
-   : face(0, 0, 0) 
+   : 
+  //face(0, 0, 0),      // NBN: replaced this MFace object
+    face(NULL),         // NBN: with a pointer to MFace
+    parentElement(0),   // NBN: also, init members
+    parentFace(0),
+    faceIndex(0),
+    zoneIndex(0)
 { }
 
 
diff --git a/Geo/MZoneBoundary.h b/Geo/MZoneBoundary.h
index 373f88e27f5d02ad5dccbc82fc9dc7af4f711c17..04cc2d6ff72cb0e4f026403eae17c4095a28b1a5 100644
--- a/Geo/MZoneBoundary.h
+++ b/Geo/MZoneBoundary.h
@@ -75,7 +75,7 @@ struct ZoneConnectivity
     int vertexIndex1;
     int vertexIndex2;
     // Constructors
-    VertexPair()
+    VertexPair() 
       : vertex(0), vertexIndex1(0), vertexIndex2(0)
     { }
     VertexPair(MVertex *const _vertex,
@@ -242,14 +242,28 @@ class MZoneBoundary
   {
     struct FaceDataB
     {
-      FaceT face;
+      // NBN: cannot use a FaceT object in FaceVector. 
+      // class FaceT has embedded std::vector objects;
+      // custom allocator for FaceVector<T> does not call ctors,
+      // but std:: dtors will try to delete _v, _si
+      // 
+      // Simple fix: use a pointer to FaceT, then build the 
+      // FaceT object once the FaceDataB structure has been 
+      // safely added to the container (push_back)
+
+    //FaceT  face;      // NBN: FaceT contains std:: containers 
+      FaceT* face;      // NBN: use FaceT* (then init in two steps)
+
       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),
+        : 
+      //face(bFMapIt->first), 
+        face(NULL),   // NBN: need to load this after insertion into container
+        parentElement(bFMapIt->second.parentElement),
         parentFace(bFMapIt->second.parentFace),
         faceIndex(bFMapIt->second.faceIndex), zoneIndex(_zoneIndex)
       { }
@@ -271,7 +285,9 @@ class MZoneBoundary
                  // fails on some compilers (earlier versions of g++?)
       // The default constructor is required by 'set_offsets()' in
       // class 'FaceAllocator'.  This is invoked by preInit() below.
-      ZoneData() { };
+      ZoneData() 
+        : vertexIndex(0), zoneIndex(0)    // NBN: init members
+      { }
       friend class CCon::FaceAllocator<ZoneData>;
     };
     CCon::FaceVector<FaceDataB> faces;
@@ -299,6 +315,27 @@ class MZoneBoundary
 
   void clear()
   {
+    // NBN: using FaceT* so need to dealloc:
+    int icount = 0;
+    GlobalBoVertexMap::iterator itEnd = globalBoVertMap.end();
+    for (GlobalBoVertexMap::iterator itBoV = globalBoVertMap.begin(); 
+          itBoV != itEnd; ++itBoV)
+    {
+      // ... clear the faces
+      typename GlobalVertexData<FaceT>& ref = itBoV->second;
+      size_t nf = ref.faces.size();
+      for (int i=0; i<nf; ++i) {
+        ++ icount;
+        FaceT* p = ref.faces[i].face;
+        if (p) {
+          delete (p);
+          p = NULL;
+        }
+      }
+    }
+    Msg::Info("cleared %d faces.", icount);
+
+    // finally, clear the container
     globalBoVertMap.clear();
   }
 
diff --git a/Plugin/MathEval.h b/Plugin/MathEval.h
index fd68e5ce77eaaeee8de3c5600f6475130af38b12..df6ff53695aa376e8b75d86415d539bab9f70aa0 100644
--- a/Plugin/MathEval.h
+++ b/Plugin/MathEval.h
@@ -7,7 +7,7 @@
 #define _MATH_EVAL_H_
 
 #include "Plugin.h"
-;
+
 extern "C"
 {
   GMSH_Plugin *GMSH_RegisterMathEvalPlugin();