diff --git a/src/mmc_highorder.cpp b/src/mmc_highorder.cpp
index b3c85ab..1d1e852 100644
--- a/src/mmc_highorder.cpp
+++ b/src/mmc_highorder.cpp
@@ -100,47 +100,33 @@ void mesh_10nodetet(tetmesh* mesh, mcconfig* cfg) {
     for (it = edgelist.begin(); it != edgelist.end(); it++) {
         for (int i = 0; i < 3; i++) {
             ((float*)(&mesh->node[oldnn + pos]))[i] =
-                (((float*)(&mesh->node[(*it).first]))[i] + ((float*)(&mesh->node[(*it).second]))[i]) * 0.5f;
+                (((float*)(&mesh->node[it->first]))[i] + ((float*)(&mesh->node[it->second]))[i]) * 0.5f;
         }
 
         pos++;
     }
 }
 
-std::string face_to_string(int a, int b, int c) {
-    int ar[3];
-    ar[0] = MIN(MIN(a, b), c);
-    ar[2] = MAX(MAX(a, b), c);
-    ar[1] = - ar[0] - ar[2] + a + b + c;
-
-    std::ostringstream strkey;
-    strkey << ar[0] << "," << ar[1] << "," << ar[2];
-
-    return std::string(strkey.str());
-}
 
 #ifdef __cplusplus
     extern "C"
 #endif
 void mesh_getfacenb(tetmesh* mesh, mcconfig* cfg) {
-    std::vector<std::string> faces;
-    std::map< std::string, std::pair<unsigned int, unsigned int> > facenb;
-
-    faces.reserve(mesh->ne << 2);
+    std::map< std::vector<int>, std::pair<unsigned int, unsigned int> > facenb;
+    std::map< std::vector<int>, std::pair<unsigned int, unsigned int> > ::iterator it;
 
     for (int i = 0; i < mesh->ne; i++) {
         int* ee = mesh->elem + i * mesh->elemlen;
 
         for (int j = 0; j < 4; j++) {
-            faces.push_back(face_to_string(ee[facelist[j][0]], ee[facelist[j][1]], ee[facelist[j][2]]));
-        }
-    }
+            std::vector<int> facevec = { ee[facelist[j][0]], ee[facelist[j][1]], ee[facelist[j][2]] };
+            std::sort (facevec.begin(), facevec.end());
 
-    for (unsigned int i = 0; i < faces.size(); i++) {
-        if (facenb.count(faces[i])) {
-            facenb[faces[i]].second = (i >> 2) + 1;
-        } else {
-            facenb[faces[i]] = std::make_pair((i >> 2) + 1, 0);
+            if (facenb.count(facevec)) {
+                facenb[facevec].second = (i << 2) + j + 1;
+            } else {
+                facenb[facevec] = std::make_pair((i << 2) + j, 0);
+            }
         }
     }
 
@@ -150,11 +136,10 @@ void mesh_getfacenb(tetmesh* mesh, mcconfig* cfg) {
 
     mesh->facenb = (int*)calloc(sizeof(int) * mesh->elemlen, mesh->ne);
 
-    for (unsigned int i = 0; i < faces.size(); i++) {
-        std::pair<unsigned int, unsigned int> item = facenb[faces[i]];
-
-        if (item.second != 0) {
-            mesh->facenb[i] = ((i >> 2) + 1 == item.first ? item.second : item.first);
+    for (it = facenb.begin(); it != facenb.end(); it++) {
+        if (it->second.second != 0) {
+            mesh->facenb[it->second.first] = ((it->second.second - 1) >> 2) + 1;
+            mesh->facenb[it->second.second - 1] = (it->second.first >> 2) + 1;
         }
     }
 }
diff --git a/src/mmc_mesh.c b/src/mmc_mesh.c
index b2e1812..295ef71 100644
--- a/src/mmc_mesh.c
+++ b/src/mmc_mesh.c
@@ -731,15 +731,11 @@ void mesh_loadelemvol(tetmesh* mesh, mcconfig* cfg) {
     int tmp, len, i, j, *ee;
     char fvelem[MAX_FULL_PATH];
 
-    if (cfg->elem && cfg->elemnum && cfg->node && cfg->nodenum) {
-        mesh_getvolume(mesh, cfg);
-        return;
-    }
-
     mesh_filenames("velem_%s.dat", fvelem, cfg);
 
-    if ((fp = fopen(fvelem, "rt")) == NULL) {
-        MESH_ERROR("can not open element volume file");
+    if ((cfg->elem && cfg->elemnum && cfg->node && cfg->nodenum) || (fp = fopen(fvelem, "rt")) == NULL) {
+        mesh_getvolume(mesh, cfg);
+        return;
     }
 
     len = fscanf(fp, "%d %d", &tmp, &(mesh->ne));
@@ -784,15 +780,11 @@ void mesh_loadfaceneighbor(tetmesh* mesh, mcconfig* cfg) {
     int* pe;
     char ffacenb[MAX_FULL_PATH];
 
-    if (cfg->elem && cfg->elemnum) {
-        mesh_getfacenb(mesh, cfg);
-        return;
-    }
-
     mesh_filenames("facenb_%s.dat", ffacenb, cfg);
 
-    if ((fp = fopen(ffacenb, "rt")) == NULL) {
-        MESH_ERROR("can not open face-neighbor list file");
+    if ((cfg->elem && cfg->elemnum) || (fp = fopen(ffacenb, "rt")) == NULL) {
+        mesh_getfacenb(mesh, cfg);
+        return;
     }
 
     len = fscanf(fp, "%d %d", &(mesh->elemlen), &(mesh->ne));
@@ -927,22 +919,23 @@ int mesh_initelem(tetmesh* mesh, mcconfig* cfg) {
     int i, j;
 
     for (i = 0; i < mesh->ne; i++) {
-        float3 pmin = {VERY_BIG, VERY_BIG, VERY_BIG}, pmax = {-VERY_BIG, -VERY_BIG, -VERY_BIG};
+        double pmin[3] = {VERY_BIG, VERY_BIG, VERY_BIG}, pmax[3] = {-VERY_BIG, -VERY_BIG, -VERY_BIG};
         int* elems = (int*)(mesh->elem + i * mesh->elemlen); // convert int4* to int*
 
         for (j = 0; j < mesh->elemlen; j++) {
-            pmin.x = MIN(nodes[elems[0] - 1].x, pmin.x);
-            pmin.y = MIN(nodes[elems[0] - 1].y, pmin.y);
-            pmin.z = MIN(nodes[elems[0] - 1].z, pmin.z);
+            pmin[0] = MIN(nodes[elems[j] - 1].x, pmin[0]);
+            pmin[1] = MIN(nodes[elems[j] - 1].y, pmin[1]);
+            pmin[2] = MIN(nodes[elems[j] - 1].z, pmin[2]);
 
-            pmax.x = MAX(nodes[elems[0] - 1].x, pmax.x);
-            pmax.y = MAX(nodes[elems[0] - 1].y, pmax.y);
-            pmax.z = MAX(nodes[elems[0] - 1].z, pmax.z);
+            pmax[0] = MAX(nodes[elems[j] - 1].x, pmax[0]);
+            pmax[1] = MAX(nodes[elems[j] - 1].y, pmax[1]);
+            pmax[2] = MAX(nodes[elems[j] - 1].z, pmax[2]);
         }
 
-        if (cfg->srcpos.x <= pmax.x && cfg->srcpos.x >= pmin.x &&
-                cfg->srcpos.y <= pmax.y && cfg->srcpos.y >= pmin.y &&
-                cfg->srcpos.z <= pmax.z && cfg->srcpos.z >= pmin.z) {
+        if (cfg->srcpos.x <= pmax[0] && cfg->srcpos.x >= pmin[0] &&
+                cfg->srcpos.y <= pmax[1] && cfg->srcpos.y >= pmin[1] &&
+                cfg->srcpos.z <= pmax[2] && cfg->srcpos.z >= pmin[2]) {
+
             if (mesh_barycentric(i + 1, &(cfg->bary0.x), (FLOAT3*) & (cfg->srcpos), mesh) == 0) {
                 cfg->e0 = i + 1;
                 return 0;