Index: XIOS/dev/branch_openmp/src/node/mesh.cpp
===================================================================
--- XIOS/dev/branch_openmp/src/node/mesh.cpp	(revision 1328)
+++ XIOS/dev/branch_openmp/src/node/mesh.cpp	(revision 1545)
@@ -7,4 +7,5 @@
 #include "mesh.hpp"
 using namespace ep_lib;
+#include <boost/functional/hash.hpp>
 
 namespace xios {
@@ -32,7 +33,4 @@
   }
 
-  //std::map <StdString, CMesh> CMesh::meshList = std::map <StdString, CMesh>();
-  //std::map <StdString, vector<int> > CMesh::domainList = std::map <StdString, vector<int> >();
-
   std::map <StdString, CMesh> *CMesh::meshList_ptr = 0;
   std::map <StdString, vector<int> > *CMesh::domainList_ptr = 0;
@@ -50,22 +48,16 @@
     if(CMesh::meshList_ptr == NULL)   CMesh::meshList_ptr   = new std::map <StdString, CMesh>();
 
-    //CMesh::domainList[meshName].push_back(nvertex);
     CMesh::domainList_ptr->at(meshName).push_back(nvertex);
 
-    //if ( CMesh::meshList.begin() != CMesh::meshList.end() )
     if ( CMesh::meshList_ptr->begin() != CMesh::meshList_ptr->end() )
     {
-      //for (std::map<StdString, CMesh>::iterator it=CMesh::meshList.begin(); it!=CMesh::meshList.end(); ++it)
       for (std::map<StdString, CMesh>::iterator it=CMesh::meshList_ptr->begin(); it!=CMesh::meshList_ptr->end(); ++it)
       {
         if (it->first == meshName)
-          //return &meshList[meshName];
           return &meshList_ptr->at(meshName);
         else
         {
           CMesh newMesh;
-          //CMesh::meshList.insert( make_pair(meshName, newMesh) );
           CMesh::meshList_ptr->insert( make_pair(meshName, newMesh) );
-          //return &meshList[meshName];
           return &meshList_ptr->at(meshName);
         }
@@ -75,7 +67,5 @@
     {
       CMesh newMesh;
-      //CMesh::meshList.insert( make_pair(meshName, newMesh) );
       CMesh::meshList_ptr->insert( make_pair(meshName, newMesh) );
-      //return &meshList[meshName];
       return &meshList_ptr->at(meshName);
     }
@@ -151,80 +141,4 @@
     }
 
-
-///----------------------------------------------------------------
-/*!
- * \fn size_t CMesh::nodeIndex (double lon, double lat)
- * Returns its index if a node exists; otherwise adds the node and returns -1.
- * Precision check is implemented with two hash values for each dimension, longitude and latitude.
- * \param [in] lon Node longitude in degrees.
- * \param [in] lat Node latitude in degrees ranged from 0 to 360.
- * \return node index if a node exists; -1 otherwise
- */
-  size_t CMesh::nodeIndex (double lon, double lat)
-  {
-    double minBoundLon = 0. ;
-    double maxBoundLon = 360. ;
-    double minBoundLat = -90 ;
-    double maxBoundLat = 90 ;
-    double prec=1e-11 ;
-    double precLon=prec ;
-    double precLat=prec ;
-
-    size_t maxsize_t=numeric_limits<size_t>::max() ;
-    if ( (maxBoundLon-minBoundLon)/maxsize_t > precLon) precLon=(maxBoundLon-minBoundLon)/maxsize_t ;
-    if ( (maxBoundLat-minBoundLat)/maxsize_t > precLat) precLat=(maxBoundLat-minBoundLat)/maxsize_t ;
-
-    size_t iMinLon=0 ;
-    size_t iMaxLon=(maxBoundLon-minBoundLon)/precLon ;
-    size_t iMinLat=0 ;
-    size_t iMaxLat=(maxBoundLat-minBoundLat)/precLat ;
-
-    size_t hash0,hash1,hash2,hash3 ;
-    size_t lon0,lon1,lat0,lat1 ;
-
-    lon0=(lon-minBoundLon)/precLon ;
-    if ( ((lon0+1)*precLon + lon0*precLon)/2 > lon-minBoundLon)
-    {
-      if (lon0==iMinLon) lon1=iMaxLon ;
-      else lon1=lon0-1 ;
-    }
-    else
-    {
-      if (lon0==iMaxLon) lon1=iMinLon ;
-      else lon1=lon0+1 ;
-    }
-
-    lat0=(lat-minBoundLat)/precLat ;
-    if ( ((lat0+1)*precLat + lat0*precLat)/2 > lat-minBoundLat)
-    {
-      if (lat0==iMinLat) lat1=lat0 ;
-      else lat1=lat0-1 ;
-    }
-    else
-    {
-      if (lat0==iMaxLat) lat1=lat0 ;
-      else lat1=lat0+1 ;
-    }
-
-    hash0=hashPair(lon0,lat0) ;
-    hash1=hashPair(lon0,lat1) ;
-    hash2=hashPair(lon1,lat0) ;
-    hash3=hashPair(lon1,lat1) ;
-
-    boost::unordered_map<size_t, size_t>::iterator end = hashed_map_nodes.end() ;
-    size_t mapSize = hashed_map_nodes.size();
-    if (hashed_map_nodes.find(hash0)==end && hashed_map_nodes.find(hash1)==end && hashed_map_nodes.find(hash2)==end && hashed_map_nodes.find(hash3)==end)
-    {
-      hashed_map_nodes[hash0] = mapSize ;
-      hashed_map_nodes[hash1] = mapSize + 1;
-      hashed_map_nodes[hash2] = mapSize + 2;
-      hashed_map_nodes[hash3] = mapSize + 3;
-      return -1;
-    }
-    else
-      return ( (hashed_map_nodes[hash0]+1) / 4 );
-
-  } // nodeIndex()
-
 ///----------------------------------------------------------------
 /*!
@@ -313,181 +227,4 @@
  * \param [in] bounds_lat Array of boundary latitudes. Its size depends on the element type.
  */
-  void CMesh::createMesh(const CArray<double, 1>& lonvalue, const CArray<double, 1>& latvalue,
-            const CArray<double, 2>& bounds_lon, const CArray<double, 2>& bounds_lat)
-  {
-    int nvertex = (bounds_lon.numElements() == 0) ? 1 : bounds_lon.rows();
-
-    if (nvertex == 1)
-    {
-      nbNodes_ = lonvalue.numElements();
-      node_lon.resizeAndPreserve(nbNodes_);
-      node_lat.resizeAndPreserve(nbNodes_);
-      for (int nn = 0; nn < nbNodes_; ++nn)
-      {
-        if (map_nodes.find(make_pair (lonvalue(nn), latvalue(nn))) == map_nodes.end())
-        {
-          map_nodes[make_pair (lonvalue(nn), latvalue(nn))] = nn ;
-          node_lon(nn) = lonvalue(nn);
-          node_lat(nn) = latvalue(nn);
-        }
-      }
-    }
-    else if (nvertex == 2)
-    {
-      nbEdges_ = bounds_lon.shape()[1];
-
-      // Create nodes and edge_node connectivity
-      node_lon.resizeAndPreserve(nbEdges_*nvertex); // Max possible number of nodes
-      node_lat.resizeAndPreserve(nbEdges_*nvertex);
-      edge_nodes.resizeAndPreserve(nvertex, nbEdges_);
-
-      for (int ne = 0; ne < nbEdges_; ++ne)
-      {
-        for (int nv = 0; nv < nvertex; ++nv)
-        {
-          if (map_nodes.find(make_pair (bounds_lon(nv, ne), bounds_lat(nv ,ne))) == map_nodes.end())
-          {
-            map_nodes[make_pair (bounds_lon(nv, ne), bounds_lat(nv, ne))] = nbNodes_ ;
-            edge_nodes(nv,ne) = nbNodes_ ;
-            node_lon(nbNodes_) = bounds_lon(nv, ne);
-            node_lat(nbNodes_) = bounds_lat(nv, ne);
-            ++nbNodes_ ;
-          }
-          else
-            edge_nodes(nv,ne) = map_nodes[make_pair (bounds_lon(nv, ne), bounds_lat(nv ,ne))];
-        }
-      }
-      node_lon.resizeAndPreserve(nbNodes_);
-      node_lat.resizeAndPreserve(nbNodes_);
-
-      // Create edges
-      edge_lon.resizeAndPreserve(nbEdges_);
-      edge_lat.resizeAndPreserve(nbEdges_);
-
-      for (int ne = 0; ne < nbEdges_; ++ne)
-      {
-        if (map_edges.find(make_ordered_pair (edge_nodes(0,ne), edge_nodes(1,ne))) == map_edges.end())
-        {
-          map_edges[make_ordered_pair ( edge_nodes(0,ne), edge_nodes(1,ne) )] = ne ;
-          edge_lon(ne) = lonvalue(ne);
-          edge_lat(ne) = latvalue(ne);
-        }
-
-      }
-      edgesAreWritten = true;
-    }
-    else
-    {
-      nbFaces_ = bounds_lon.shape()[1];
-  
-      // Create nodes and face_node connectivity
-      node_lon.resizeAndPreserve(nbFaces_*nvertex);  // Max possible number of nodes
-      node_lat.resizeAndPreserve(nbFaces_*nvertex);
-      face_nodes.resize(nvertex, nbFaces_);
-  
-      for (int nf = 0; nf < nbFaces_; ++nf)
-      {
-        for (int nv = 0; nv < nvertex; ++nv)
-        {
-          if (map_nodes.find(make_pair (bounds_lon(nv, nf), bounds_lat(nv ,nf))) == map_nodes.end())
-          {
-            map_nodes[make_pair (bounds_lon(nv, nf), bounds_lat(nv, nf))] = nbNodes_ ;
-            face_nodes(nv,nf) = nbNodes_ ;
-            node_lon(nbNodes_) = bounds_lon(nv, nf);
-            node_lat(nbNodes_) = bounds_lat(nv ,nf);
-            ++nbNodes_ ;
-          }
-          else
-          {
-            face_nodes(nv,nf) = map_nodes[make_pair (bounds_lon(nv, nf), bounds_lat(nv ,nf))];
-          }
-        }
-      }
-      node_lon.resizeAndPreserve(nbNodes_);
-      node_lat.resizeAndPreserve(nbNodes_);
-  
-      // Create edges and edge_nodes connectivity
-      edge_lon.resizeAndPreserve(nbFaces_*nvertex); // Max possible number of edges
-      edge_lat.resizeAndPreserve(nbFaces_*nvertex);
-      edge_nodes.resizeAndPreserve(2, nbFaces_*nvertex);
-      edge_faces.resize(2, nbFaces_*nvertex);
-      face_edges.resize(nvertex, nbFaces_);
-      face_faces.resize(nvertex, nbFaces_);
-
-      vector<int> countEdges(nbFaces_*nvertex);   // needed in case if edges have been already generated
-      vector<int> countFaces(nbFaces_);
-      countEdges.assign(nbFaces_*nvertex, 0);
-      countFaces.assign(nbFaces_, 0);
-      int edge;
-      for (int nf = 0; nf < nbFaces_; ++nf)
-      {
-        for (int nv1 = 0; nv1 < nvertex; ++nv1)
-        {
-          int nv = 0;
-          int nv2 = (nv1 < nvertex -1 ) ? (nv1 + 1) : (nv1 + 1 - nvertex); // cyclic rotation
-          if (map_edges.find(make_ordered_pair (face_nodes(nv1,nf), face_nodes(nv2,nf))) == map_edges.end())
-          {
-            map_edges[make_ordered_pair (face_nodes(nv1,nf), face_nodes(nv2,nf))] = nbEdges_ ;
-            face_edges(nv1,nf) = map_edges[make_ordered_pair (face_nodes(nv1,nf), face_nodes(nv2,nf))];
-            edge_faces(0,nbEdges_) = nf;
-            edge_faces(1,nbEdges_) = -999;
-            face_faces(nv1,nf) = 999999;
-            edge_nodes(Range::all(),nbEdges_) = face_nodes(nv1,nf), face_nodes(nv2,nf);
-            edge_lon(nbEdges_) = ( abs( node_lon(face_nodes(nv1,nf)) - node_lon(face_nodes(nv2,nf))) < 180.) ?
-                        (( node_lon(face_nodes(nv1,nf)) + node_lon(face_nodes(nv2,nf))) * 0.5) :
-                        (( node_lon(face_nodes(nv1,nf)) + node_lon(face_nodes(nv2,nf))) * 0.5 -180.);
-            edge_lat(nbEdges_) = ( node_lat(face_nodes(nv1,nf)) + node_lat(face_nodes(nv2,nf)) ) * 0.5;
-            ++nbEdges_;
-          }
-          else
-          {
-            edge = map_edges[make_ordered_pair (face_nodes(nv1,nf), face_nodes(nv2,nf))];
-            face_edges(nv1,nf) = edge;
-            if (edgesAreWritten)
-            {
-              edge_faces(countEdges[edge], edge) = nf;
-              if (countEdges[edge]==0)
-              {
-                face_faces(nv1,nf) = 999999;
-              }
-              else
-              {
-                int face1 = nf; // = edge_faces(1,edge)
-                int face2 = edge_faces(0,edge);
-                face_faces(countFaces[face1], face1) =  face2;
-                face_faces(countFaces[face2], face2) =  face1;
-                ++(countFaces[face1]);
-                ++(countFaces[face2]);
-              }
-            }
-            else
-            {
-              edge_faces(1,edge) = nf;
-              int face1 = nf; // = edge_faces(1,edge)
-              int face2 = edge_faces(0,edge);
-              face_faces(countFaces[face1], face1) =  face2;
-              face_faces(countFaces[face2], face2) =  face1;
-              ++(countFaces[face1]);
-              ++(countFaces[face2]);
-            }
-            ++(countEdges[edge]);
-          }
-        }
-      }
-      edge_nodes.resizeAndPreserve(2, nbEdges_);
-      edge_faces.resizeAndPreserve(2, nbEdges_);
-      edge_lon.resizeAndPreserve(nbEdges_);
-      edge_lat.resizeAndPreserve(nbEdges_);
-
-      // Create faces
-      face_lon.resize(nbFaces_);
-      face_lat.resize(nbFaces_);
-      face_lon = lonvalue;
-      face_lat = latvalue;
-      facesAreWritten = true;
-
-    } // nvertex > 2
-    
-  } // createMesh()
 
 ///----------------------------------------------------------------
@@ -547,6 +284,6 @@
 
       // For determining the global edge index
-      size_t nbEdgesOnProc = nbEdges_;
-      size_t nbEdgesAccum;
+      unsigned long nbEdgesOnProc = nbEdges_;
+      unsigned long nbEdgesAccum;
       MPI_Scan(&nbEdgesOnProc, &nbEdgesAccum, 1, MPI_UNSIGNED_LONG, MPI_SUM, comm);
       nbEdgesAccum -= nbEdges_;
@@ -677,6 +414,6 @@
          // Maps modified in this step:
          // nodeIdx2Idx = <idx, idxGlo>
-         int nodeCount = nodeIdx2Idx.size();
-         int nodeStart, nbNodes;
+         unsigned long nodeCount = nodeIdx2Idx.size();
+         unsigned long nodeStart, nbNodes;
          MPI_Scan(&nodeCount, &nodeStart, 1, MPI_UNSIGNED_LONG, MPI_SUM, comm);
          int nNodes = nodeStart;
@@ -770,6 +507,6 @@
 
       // For determining the global face index
-      size_t nbFacesOnProc = nbFaces_;
-      size_t nbFacesAccum;
+      unsigned long nbFacesOnProc = nbFaces_;
+      unsigned long nbFacesAccum;
       MPI_Scan(&nbFacesOnProc, &nbFacesAccum, 1, MPI_UNSIGNED_LONG, MPI_SUM, comm);
       nbFacesAccum -= nbFaces_;
@@ -895,5 +632,5 @@
         }
 
-        int edgeStart, nbEdges;
+        unsigned long edgeStart, nbEdges;
         MPI_Scan(&edgeCount, &edgeStart, 1, MPI_UNSIGNED_LONG, MPI_SUM, comm);
         int nEdges = edgeStart;
@@ -1115,6 +852,6 @@
         }
 
-        int edgeCount = edgeIdx2Idx.size();
-        int edgeStart, nbEdges;
+        unsigned long edgeCount = edgeIdx2Idx.size();
+        unsigned long edgeStart, nbEdges;
         MPI_Scan(&edgeCount, &edgeStart, 1, MPI_UNSIGNED_LONG, MPI_SUM, comm);
         int nEdges = edgeStart;
@@ -1385,6 +1122,6 @@
         // Maps modified in this step:
         // nodeIdx2Idx = <idx, idxGlo>
-        int nodeCount = nodeIdx2Idx.size();
-        int nodeStart, nbNodes;
+        unsigned long nodeCount = nodeIdx2Idx.size();
+        unsigned long nodeStart, nbNodes;
         MPI_Scan(&nodeCount, &nodeStart, 1, MPI_UNSIGNED_LONG, MPI_SUM, comm);
         int nNodes = nodeStart;
@@ -1505,6 +1242,6 @@
         }
 
-        int edgeCount = edgeIdx2Idx.size();
-        int edgeStart, nbEdges;
+        unsigned long edgeCount = edgeIdx2Idx.size();
+        unsigned long edgeStart, nbEdges;
         MPI_Scan(&edgeCount, &edgeStart, 1, MPI_UNSIGNED_LONG, MPI_SUM, comm);
         int nEdges = edgeStart;
@@ -2103,5 +1840,5 @@
 
     // faceToFaces connectivity
-    boost::unordered_map <int, int> mapFaces;  // mapFaces = < hash(face1, face2), hash> (the mapped value is irrelevant)
+    std::unordered_map <int, int> mapFaces;  // mapFaces = < hash(face1, face2), hash> (the mapped value is irrelevant)
     int maxNb = 20;                            // some assumption on the max possible number of neighboring cells
     faceToFaces.resize(maxNb, nbFaces);
@@ -2160,5 +1897,5 @@
     CArray<double, 2> faceToNodes (nvertex, nbFaces);
 
-    boost::unordered_map <pair<double,double>, int> mapNodes;
+    std::unordered_map <pairDouble, int, boost::hash<pairDouble> > mapNodes;
 
     for (int nf = 0; nf < nbFaces; ++nf)
@@ -2176,5 +1913,5 @@
 
     // faceToFaces connectivity
-    boost::unordered_map <pair<int,int>, int> mapEdges;
+    std::unordered_map <pairInt, int, boost::hash<pairInt> > mapEdges;
     faceToFaces.resize(nvertex, nbFaces);
     CArray<int, 2> edgeToFaces(2, nbFaces*nvertex); // max possible
